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1  Introduction 

1.1  Background 

Breast  cancer  is  the  most  frequently  diagnosed  invasive  malignancy  among 
American  women  and  ranks  second  only  to  lung  cancer  in  annual  cancer-related 
mortality  for  this  group  [1].  Numerous  studies  have  shown  that  early  detection 
and  treatment  of  breast  cancer  can  improve  survival  rates  [2-4].  Screen-film 
mammography  has  come  to  play  a  vital  role  in  this  detection  process,  due  to  its 
high  (80-90%)  sensitivity  to  breast  malignancies.  However,  mammography  is  no¬ 
toriously  poor  at  distinguishing  benign  from  malignant  tumors,  having  reported 
specificities  and  positive  predictive  values  of  15-30%  [5].  This  means  that  only 
a  small  percentage  of  lesions  biopsied  on  the  basis  of  suspicious  mammographic 
appearance  are  found  to  be  malignant. 

In  recent  years,  researchers  have  developed  and  studied  a  nuclear-medicine 
test  with  the  potential  to  provide  relatively  low-cost,  minimally  invasive  differ¬ 
entiation  of  breast  abnormalities  identified  by  physical  examination  or  mam¬ 
mography  [6-18].  Known  as  scintimammography  (SMM),  the  test  relies  on  the 
preferential  uptake  of  Tc-99m-sestamibi  or  other  radionuclides  such  as  Tl-201, 
Tc-99m-tetrofosmin,  or  Tc-99m-MDP  in  breast  malignancies  as  compared  to 
normal  breast  tissue  or  benign  abnormalities.  Indeed,  one  study  has  shown  that 
typical  in  vivo  tumor-background  concentration  ratios  of  Tc-99m-sestamibi  are 
5.64±3.06  [19].  This  focal  uptake  can  be  imaged  in  a  number  of  ways,  though 
the  most  widely  used  clinical  protocol  involves  acquiring  one  or  two  planar 
views — one  lateral  view  and  possibly  an  additional  oblique  or  anterior  view — 
while  the  patient  lies  prone  on  a  specially  designed  table  [20].  The  imaging  time 
is  typically  10-15  minutes  per  view.  Numerous  clinical  studies  with  histological 
follow-up  have  been  performed  using  this  or  a  similar  protocol,  reporting  sensi¬ 
tivities  of  83-96%,  and  specificities  of  66-100%  when  using  Tc-99m-sestamibi  to 
image  mammographically  suspicious  lesions  [8-18].  In  addition  to  differentiat¬ 
ing  breast  abnormalities  detected  by  other  means,  SMM  may  also  have  a  role 
in  the  detection  of  breast  malignancies  in  patients  with  radiographically  dense 
breasts,  for  whom  screen-film  mammograms  are  often  difficult  to  interpret  [7]. 

A  few  of  these  studies  of  SMM  have  also  examined  the  role  of  conventional 
single  photon  emission  computed  tomography  (SPECT)  (where  the  patient  lies 
supine  and  the  camera  revolves  around  the  torso)  in  detecting  focal  uptake 
of  Tc-99m-sestamibi  and  have  found  comparable  but  not  generally  improved 
sensitivities  as  compared  to  planar  techniques,  along  with  substantially  lower 
specificities  [14,21-23].  Wang  et  al  [24]  speculated  that  this  surprisingly  poor 
performance  was  due  to  substantial  attenuation  of  photons  emitted  in  the  breast 
by  the  torso  in  at  least  half  of  the  views  as  well  as  to  the  presence  of  scatter 
from  organs  such  as  the  heart  and  liver  known  to  have  high  uptake  of  Tc-99m- 
sestamibi.  The  poor  performance  of  conventional  SPECT  may  also  be  related  to 
the  inferior  resolution  of  conventional  SPECT  as  compared  to  planar  techniques 
in  this  situation,  due  to  the  fact  that  the  scintillation  cameras  are  on  average 
further  away  from  the  breast  in  the  conventional  SPECT  geometry  than  in  the 
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planar  geometry.  Wang  et  al  also  investigated  a  geometry  that  they  called 
vertical  axis-of-rotation  SPECT  and  that  we  call  dedicated  breast  SPECT,  or 
simply  dedicated  SPECT,  in  which  the  scintillation  cameras  are  assumed  to 
revolve  around  one  breast  alone,  with  the  patient  lying  prone.  This  geometry 
eliminates  the  effect  of  attenuation  by  the  thorax  and,  with  proper  shielding,  the 
effect  of  scatter  from  the  thorax.  Moreover,  the  small  radius  of  rotation  offers 
improved  resolution  and  sensitivity.  In  phantom  studies,  Wang  et  al  found 
that  with  this  dedicated  geometry  they  were  able  to  detect  a  breast  lesion  with 
an  outer  diameter  of  1  cm  and  a  6:1  lesion-to-background  concentration  ratio 
that  was  not  detectable  in  either  conventional  SPECT  or  planar  studies  with 
the  same  total  imaging  time. 

While  no  one  has  yet  developed  a  dedicated  breast  SPECT  system  for  clin¬ 
ical  use,  encouraging  results  such  as  these  make  it  an  important  and  hopefully 
fruitful  area  of  research.  One  practical  concern  that  could  limit  the  acceptance 
and  use  of  dedicated  SPECT  SMM  is  the  lengthy  imaging  time  (on  the  order 
of  30-40  minutes)  that  would  be  required  to  acquire  projection  data  useful  to 
standard  reconstruction  algorithms  such  as  filtered  backprojection  (EBP).  Mo¬ 
tion  artifacts  could  be  common  as  patients  could  find  it  difficult  to  remain  still 
for  this  time  while  lying  prone  on  an  unpadded  table  (padding  would  interfere 
with  imaging  near  the  chest  wall).  Moreover,  this  lengthy  imaging  time  would 
limit  patient  throughput,  which  is  undesirable  for  a  test  that  could  potentially 
gain  high- volume  use  as  a  follow-up  to  conventional  mammography. 

1.2  Purpose  of  research 

The  purpose  of  the  research  being  undertaken  is  to  reduce  imaging  time  for  ded¬ 
icated  SPECT  scintimammography  by  developing  algorithms  tailored  to  gener¬ 
ate  diagnostically  useful  images  from  a  smaller  number  of  projection  views  than 
are  usually  used  in  emission  tomography.  The  algorithms  will  take  advantage 
of  the  simplicities  and  symmetries  inherent  in  the  SMM  problem  to  extract  a 
maximum  of  information  from  the  views  that  are  acquired,  in  contrast  to  recon¬ 
struction  approaches  such  as  the  ever-popular  filtered  backprojection  (EBP), 
which  generally  force  the  data  to  yield  to  their  own  implicit  assumptions  and 
lead  to  artifacts  when  they  are  not  satisfied.  Eor  example,  it  is  well  known  that 
EBP  yields  images  with  prominent  star  artifacts  when  the  number  of  angular 
views  is  insufficient. 

In  the  proposal,  we  proposed  to  investigate  the  use  of  a  class  of  reconstruction 
approaches  known  as  Algebraic  Reconstruction  Technique  (ART)  algorithms, 
which  make  no  implicit  assumptions  about  the  completeness  of  the  dataset  but 
rather  iteratively  update  the  image  in  search  of  the  solution  that  best  agrees 
with  the  data  (the  precise  definition  of  “best”  allows  for  many  variants  of  the 
basic  ART  algorithm).  Specifically,  the  work  of  the  first  year  was  to  have  met 
three  principal  goals,  as  outlined  in  the  proposal’s  Statement  of  Work: 

1.  To  finish  implementing  the  standard  ART  algorithm. 
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2,  To  systematically  explore  the  dependence  of  reconstruction  image  quality 
on  the  number  of  projections  and  the  complexity  of  the  input  image. 

3.  To  modify  the  ART  algorithm  to  perform  a  purely  binary  reconstruction 
of  the  tracer  distribution. 

1.3  Overview  of  research  accomplishments 

In  practice,  the  year’s  research  began  by  examining  a  more  basic  issue  than 
that  of  optimizing  reconstruction  algorithms:  the  question  of  whether  a  dedi¬ 
cated  SPECT  geometry  would  indeed  be  better  for  SMM  lesion  detection  than 
currently  existing  planar  or  conventional  SPECT  geometries.  The  experimental 
results  of  Wang  et  ah  [24]  mentioned  above  suggested  that  it  would  be,  but 
we  undertook  to  answer  the  question  more  definitively  and  quantitatively  by 
using  the  so-called  ideal-observer  framework  to  quantify  the  amount  of  informa¬ 
tion  contained  in  the  projections  of  a  breast  phantom  using  the  three  different 
geometries.  The  results  indicated  conclusively  that  the  dedicated  SPECT  ge¬ 
ometry  is  superior.  The  same  framework  was  then  used  to  address  goal  2 — the 
determination  of  the  dependence  of  reconstructed  image  quality  on  the  number 
of  projections  acquired. 

After  firmly  establishing  the  superiority  of  the  dedicated  SPECT  geometry, 
we  returned  to  the  question  of  developing  reconstruction  algorithms  tailored  to 
few- view  SMM  data.  Goal  1 — implementation  of  the  standard  ART  algorithm — 
was  accomplished,  though  with  disappointing  results.  Even  in  the  absence  of 
noise,  the  quality  of  the  reconstructed  images  as  compared  to  EBP  was  not 
sufficiently  high  to  justify  the  large  computational  burden  of  the  approach.  In 
the  presence  of  noise,  images  reconstructed  by  the  ART  approach  were  really 
quite  poor,  even  compared  to  ramp-filtered  FBR  These  same  general  trends 
held  for  reconstructions  from  standard  and  small  numbers  of  views,  though 
the  few-view  EBP  reconstructions  were  of  course  plagued  by  the  star  artifacts 
mentioned  in  Section  1.2.  This  poor  performance  of  ART  suggested  that  the 
exploration  of  binary  ART  reconstructions — goal  3 — would  not  be  fruitful,  and 
that  incorporation  of  the  more  explicitly  statistical  techniques  and  statistical 
priors  planned  for  the  second  year’s  work  would  be  more  appropriate.  In  this 
vein,  we  have  made  an  extensive  search  of  the  literature  on  statistical  methods 
in  image  reconstruction,  gathering  and  studying  nearly  50  articles  on  the  topic. 

The  surprisingly  strong  performance  of  EBP  that  came  to  light  in  the  course 
of  testing  the  ART  approach,  even  in  the  face  of  a  small  number  of  projections, 
suggested  that  it  might  be  fruitful  to  explore  sinogram  preprocessing  approaches 
aimed  at  reducing  the  visual  distraction  of  the  star  artifacts  while  preserving 
crucial  image  information.  Put  differently,  we  realized  that  the  relative  simplic¬ 
ity  and  symmetry  of  the  object  being  imaged  in  dedicated  SPECT  SMM  means 
that  the  number  of  angular  views  dictated  by  the  Nyquist  sampling  condition 
is  likely  to  be  smaller  than  the  number  required  to  avoid  the  severe  star  arti¬ 
facts  common  to  few-view  EBP  reconstructions.  In  this  situation,  the  additional 
views  needed  by  EBP  can  in  principal  be  interpolated  from  the  measured  ones. 
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This  led  to  a  detailed  investigation  of  the  properties  of  four  popular  inter¬ 
polation  approaches  for  periodic  functions:  linear  interpolation  (with  periodic 
boundary  conditions),  periodic  spline  interpolation,  circular  sampling  theorem 
interpolation,  and  zero-padding  interpolation.  The  investigation  consisted  of 
three  principal  elements.  First,  we  derived  previously  unappreciated  equiva¬ 
lences  among  the  methods,  showing  that  zero-padding  interpolation  and  circular 
sampling  theorem  interpolation  are  in  some  sense  mathematically  equivalent  and 
that  both  are  exact  and  equivalent  to  Whittaker-Shannon  interpolation  when 
interpolating  a  periodic,  bandlimited  function  sampled  in  accordance  wih  the 
Nyquist  sampling  condition.  Second,  we  performed  an  empirical  study  of  the 
accuracy  of  the  four  methods  in  interpolating  typical  tomographic  sinograms, 
finding  periodic  spline  interpolation  to  be  generally  superior  in  this  situation. 
Finally,  we  explored  the  noise  properties  of  the  four  methods,  deriving  analytic 
expressions  for  the  covariance  of  interpolated  points  given  that  the  measured 
samples  are  corrupted  by  additive,  zero-mean  noise  and  we  confirmed  these 
findings  with  Monte  Carlo  studies. 

Of  course,  interpolating  among  noisy  data  can  be  unwise,  regardless  of  the 
noise  properties  of  the  interpolation  method,  so  the  next  element  of  the  year’s 
work  involved  investigating  the  use  of  a  novel,  effectively  two-dimensional,  adap¬ 
tive  smoothing  algorithm  for  smoothing  the  few- view  sinogram  prior  to  interpo¬ 
lation.  The  technique  exploits  Fourier  transforms  to  reduce  the  dimensionality 
of  the  smoothing  problem  and  then  uses  the  popular  penalized  least-squares 
smoothing  measure  where  the  smoothing  parameter  is  chosen  automatically  by 
GCV.  We  studied  the  use  of  the  smoothing  technique  alone,  on  sinograms  with 
a  standard  number  of  projections,  and  also  in  conjunction  with  interpolation 
for  few-view  sinograms. 

Finally,  because  most  of  the  sinogram  preprocessing  approaches  described 
so  far  have  involved  fitting  of  continuous  smoothing  and  interpolating  splines  to 
the  measured  data,  we  explored  an  alternative  reconstruction  approach  to  FBP 
that  allows  reconstruction  to  proceed  directly  from  the  coefficients  of  the  splines 
fit  at  each  projection  angle.  The  approach  was  first  proposed  by  Wahba  [25],  but 
we  have  simplified  her  approach,  using  relationships  among  spline  coefficients  to 
eliminate  numerically  unstable  terms,  and,  perhaps  most  significantly,  extended 
the  approach  to  the  three-dimensional  Radon  inversion  problem,  where  it  takes 
on  a  much  simpler  form  than  in  two  dimensions. 

2  Body 

2.1  Comparison  of  imaging  geometries 

2.1.1  Background 

We  examined  quantitatively  the  question  of  lesion  detectability  in  three  different 
scintimammographic  geometries — planar,  conventional  SPECT,  and  dedicated 
SPECT — using  the  so-called  ideal-observer  framework  to  calculate  signal-to- 
noise  ratios  as  a  function  of  lesion  concentration  [26].  The  ideal-observer  frame- 
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work  [27]  offers  a  way  of  assessing  the  amount  of  information  the  data  from  an 
imaging  device  contain  with  regard  to  the  performance  of  a  specified  task.  For 
example,  the  simplest  such  task  is  the  detection  of  a  signal  of  known  strength, 
shape,  and  location  in  a  specified  background.  In  this  case,  the  framework  seeks 
to  quantify  the  degree  to  which  an  ideal  observer — one  who  can  use  the  informa¬ 
tion  contained  in  the  images  to  its  fullest  extent — can  reliably  distinguish  images 
containing  the  background  alone  from  images  containing  the  background  and 
the  signal  when  both  kinds  of  images  are  corrupted  by  noise,  blurring,  and  other 
imperfections.  For  linear  imaging  processes  in  which  the  noise  in  the  output  im¬ 
age  is  assumed  to  be  additive,  Gaussian,  zero-mean,  stationary,  and  independent 
of  the  presence  or  absence  of  the  signal,  the  ideal-observer  framework  allows  us 
to  characterize  fully  the  quality  of  the  imaging  system  data  with  respect  to 
the  performance  of  the  specified  signal-detection  task  by  a  single  number,  the 
ideal-observer  signal-to-noise  ratio  (SNR).  This  is  usually  expressed  as 


SNR^j 


W{v) 


dvy 


(1) 


where  K  is  the  large-scale  transfer  characteristic  of  the  imaging  system  at  the 
desired  operating  point,  MTF{v)  is  the  system  modulation  transfer  function, 
W{v)  is  the  system  Wiener  spectrum,  and  |AStn(^)|^  is  the  power  spectrum 
of  the  signal  in  input  space,  i.e.  before  it  has  been  scaled  and  degraded  by 
the  imaging  system.  This  expression  allows  one  to  determine  the  ideal-observer 
SNR  for  any  analytically  specified  input  signal  once  K,  MTF{v),  and  W{v) 
are  known.  Alternatively,  if  one  wishes  to  determine  the  ideal-observer  SNR  for 
a  particular  real  signal,  equation  (1)  could  be  re-expressed  as 


SNRj 


j 


W{v) 


(2) 


where  |A5ou«(t’)|^  is  the  power  spectrum  of  the  signal  in  output  space,  i.e.  after 
it  has  been  scaled  and  degraded  by  the  imaging  system. 

The  images  whose  information  content  we  wish  to  measure  are  not  recon¬ 
structed  images  but  rather  the  raw  projection  images  acquired  by  the  various 
imaging  systems,  for  these  offer  the  purest  measure  of  the  quality  of  the  data 
acquired  by  the  imaging  system.  The  process  of  image  reconstruction,  while 
certainly  helpful  to  human  observers,  can  never  increase  the  ideal-observer  SNR 
and  generally  reduces  it.  Using  reconstructed  images  to  compare  imaging  ge¬ 
ometries  would  thus  conflate  issues  of  data  quality  with  issues  of  reconstruction 
algorithm  accuracy.  (Of  course,  once  the  ideal-observer  SNR  of  the  projection 
data  is  known,  ideal-observer  analysis  of  reconstructed  images  can  be  used  to 
evaluate  which  reconstruction  algorithms  best  preserve  the  information  content 
of  the  raw  projections).  Because  the  noise  in  the  projection  images  is  uncorre¬ 
lated,  equation  (2)  takes  on  a  particularly  simple  form  in  this  case, 

SNR'j  =  (3) 
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where  As  is  the  signal  vector,  diag  {jrepresents  a  diagonal  matrix  and  </>j) 
is  the  noise  free  projection  of  the  background  at  projection  bin  and  projection 
angles  <f)j  [28]. 

2.1.2  Methods 

As  discussed  above,  computing  the  ideal-observer  SNR  for  a  given  imaging  ge¬ 
ometry  and  processing  approach  requires  two  ensembles  of  projection  images: 
one  set  consisting  of  images  of  the  signal  alone  and  one  set  consisting  of  images 
of  the  background  alone.  Our  background  consisted  of  an  14-cm  diameter,  800 
cc  cylindrical  phantom  and  our  signal  was  a  1-cm  outer  diameter  spherical  le¬ 
sion  insert.  This  lesion  size  is  representative  of  the  smallest  currently  detected 
in  scintimammography.  For  each  geometry  we  acquired  20  images  of  the  back¬ 
ground  alone,  for  which  we  put  3.7  mCi  of  activity  in  the  phantom  and  imaged 
for  1  minute  total  for  each  conventional  or  dedicated  SPECT  acquisition  and 
30  seconds  for  each  planar  view.  These  combinations  of  activity  and  imaging 
time  were  chosen  to  produce  clinically  realistic  count  levels  using  the  following 
reasoning.  As  with  Wang  et  al^  we  began  with  the  assumption  that  1%  of  a 
typical  25  mCi  clinical  dose  of  Tc-99m-sestamibi  is  taken  up  by  the  myocardium. 
Using  the  volume  of  the  myocardium  in  the  Data  Spectrum  Corporation  car¬ 
diac  insert  as  a  guide,  along  with  the  assumption  that  soft  tissue  will  have  a 
1:15  concentration  relative  to  the  myocardium  allowed  us  to  determine  the  ex¬ 
pected  concentration  of  activity  in  healthy  breast  tissue.  We  wished  to  compare 
detectability  in  these  three  geometries  given  the  same  total  imaging  time,  and 
assumed  that  typical  clinical  imaging  times  would  be  30  minutes  per  SPECT 
study  and  15  minutes  per  planar  view.  Thus  we  were  comparing  the  three  ge¬ 
ometries  for  equal  total  imaging  times  given  that  two-view  planar  studies  are 
common.  Given  this,  we  scaled  up  the  calculated  concentration  by  a  factor  of 
30  and  scaled  down  the  imaging  times  by  the  same  factor  to  allow  for  more 
rapid  data  acquisition.  All  imaging  times  were  adjusted  to  compensate  for  the 
decay  of  the  activity.  To  image  the  lesion  we  filled  it  with  7.6  mCi  of  Tc-99m, 
placed  it  in  the  cylinder  now  filled  with  cold  (zero-activity)  water,  and  imaged 
for  30  minutes  in  conventional  and  dedicated  SPECT  and  20  minutes  for  the 
planar  view.  This  combination  of  activity  and  imaging  time  was  chosen  simply 
to  provide  high-count,  low  noise  images  of  the  signal,  as  required  by  the  theory. 

The  dedicated  SPECT  images  were  acquired  by  placing  the  phantom  at  the 
center  of  rotation  of  a  Picker  XP2000  two-headed  SPECT  system  with  the  heads 
rotating  at  the  minimum  radius  of  rotation  (9.0  cm).  In  this  configuration  the 
heads  were  within  2.0  cm  of  the  walls  of  the  phantom.  The  breast  phantom 
was  not  attached  to  an  anthropomorphic  torso  phantom  because  Wang  et  al. 
showed  that  with  proper  shielding  the  contribution  of  scatter  from  the  torso  can 
be  made  negligible.  We  acquired  120  views  over  360°  with  each  head  acquiring 
to  a  128x128  matrix  (pixel  size=4.67  mm).  We  used  a  low-energy,  ultra-high  res¬ 
olution  collimator.  The  conventional  SPECT  images  were  also  acquired  in  the 
absence  of  an  anthropomorphic  torso  phantom,  although  the  radius  of  rotation 
(25  cm)  and  the  placement  of  the  breast  phantom  (17  cm  off-center)  were  deter- 
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mined  with  the  torso  phantom  in  place.  The  reason  for  this  curious  arrangement 
was  to  isolate  the  effect  of  the  large  radius  of  rotation  on  lesion  detectability, 
without  the  additional  degradations  caused  by  attenuation  or  scatter  in  the 
torso.  The  ideal-observer  SNR  results  for  this  arrangement  will  thus  represent 
an  upper  bound  on  the  true  detectability.  In  other  respects,  the  conventional 
SPECT  images  used  the  same  acquisition  parameters  as  the  dedicated  SPECT. 
Finally,  the  planar  view  was  acquired  with  the  phantom  flush  against  one  head, 
which  acquired  on  a  128x128  matrix  with  a  magnification  factor  of  2.0  (pixel 
size”2.33  mm).  The  ideal-observer  SNR  for  the  planar  view  was  calculated  us¬ 
ing  equation  (3)  directly.  For  the  tomographic  studies,  equation  (3)  was  applied 
to  the  two-dimensional  sinogram  corrsponding  to  the  slice  through  the  center 
of  the  lesion. 

2.1.3  Results  and  conclusions 

The  ideal-observer  SNRs  for  the  detection  of  a  1-cm  lesion  with  a  6:1  lesion- 
background  concentration  ratio  are  listed  in  Table  1  for  the  three  geometries. 
The  SNR  values  given  suggest  that  a  dedicated  SPECT  geometry  would  lead  to 
improved  detectability  for  clinically  typical  lesions  over  the  planar  and  conven¬ 
tional  SPECT  geometries,  especially  since  in  the  presence  of  attenuation  and 
scatter  from  the  torso  we  would  expect  the  difference  between  the  dedicated  and 
conventional  SPECT  geometries  to  be  even  greater  than  it  is  here.  The  success 
of  the  dedicated  geometry  can  be  attributed  to  the  fact  that  it  combines  the 
advantages  of  the  other  two  approaches:  because  of  its  small  radius  of  rotation, 
it  offers  good  sensitivity  and  resolution  comparable  to  that  of  a  planar  view 
acquired  with  the  scintillation  camera  flush  against  the  phantom,  while  it  of¬ 
fers  the  improved  contrast  offered  by  a  tomographic  system’s  ability  to  separate 
lesion  activity  from  overlying  and  underlying  activity. 


Planar  SNR 

Conventional  SPECT  SNR 

Dedicated  SPECT  SNR 

6.2 

6.7  10.7 

Table  1:  Ideal-observer  SNRs  for  detection  of  a  6:1,  1-cm  lesion  using  planar, 
conventional  SPECT,  and  dedicated  SPECT  scintimammography. 

Having  used  the  ideal-observer  SNR  framework  to  confirm  that  a  dedicated 
SPECT  geometry  should  offer  superior  lesion  detectability  to  the  conventional 
SPECT  or  planar  geometries,  we  used  the  same  approach  to  address  goal  2, 
the  dependence  of  the  dedicated  SPECT  projection  data  SNR  on  the  number 
of  input  angles  used.  The  result  is  illustrated  in  Fig,  1.  The  SNR  is  seen 
to  increase  approximately  as  the  square  root  of  the  number  of  angles  used, 
indicating  a  diminishing  SNR  return  as  the  number  of  angles  is  augmented,  and 
holding  some  promise  for  the  potential  of  few- view  imaging  techniques. 
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Figure  1:  Dependence  of  sinogram  ideal-observer  SNR  on  the  number  of  angular 
views  used  in  imaging  a  1-cm,  6:1  lesion  with  a  dedicated  SPECT  geometry.  For 
reference,  the  curve  \/Angles  is  plotted  as  well. 


2.2  Algebraic  reconstruction  techniques 

2.2.1  Background 


Having  established  the  superiority  of  the  dedicated  SPECT  geometry  and  the 
promising  square-root  dependence  of  SNR  on  the  number  of  angles  used,  we  be¬ 
gan  to  concentrate  on  the  development  of  algorithms  capable  of  producing  clini¬ 
cally  useful  images  from  a  smaller  number  of  angular  views  than  is  usually  used. 
We  focused  first  on  the  class  of  algorithms  known  as  algebraic  reconstruction 
techniques  (ART),  for  unlike  FBP  these  make  no  implicit  assumptions  about 
the  completeness  of  the  data  but  rather  seek  to  produce  iteratively  a  solution 
that  best  matches  the  data  is  available,  according  to  some  precise  definition  of 
“best.” 

In  their  most  basic  form,  ART  algorithms  begin  with  an  initial  “guess”  of  the 
image,  often  just  a  uniform  image,  then  project  this  estimate  and  compare  to 
the  measured  data.  The  resulting  differences  are  used  to  update  the  image  pixel 
values.  The  process  continues  either  for  a  prespecified  number  of  iterations  or 
until  the  magnitude  of  the  changes  from  one  iteration  to  the  next  falls  below 
some  prespecified  threshold.  One  popular  form  of  the  iterative  update  procedure 
is  given  by  the  following  equation 


U  5 


(4) 


where  and  represent  the  current  estimate  of  the  image  vector  after 

the  kth  and  (fe  H-  l)st  iterations,  respectively,  ik  =  [A:(mod/)  -h  1],  /  is  the  total 
number  of  projection  rays,  yi^is  the  measurement  for  the  ikth  projection  ray, 
is  the  transpose  of  the  uth  row  of  the  projection  matrix  R  relating  the  image 
vector  to  the  measurement  vector,  and  ()  denotes  an  inner  product  [29,  chapter 
10], 
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2.2.2  Methods 


In  accordance  with  goal  1,  mentioned  in  Section  1.2,  we  finished  implementing 
an  ART  algorithm  that  made  use  of  equation  (4).  The  iyjth  element  of  the 
projection  matrix  R  was  chosen  to  equal  the  intersection  length  between  the 
ith  pixel  and  the  jth  projection  ray.  To  calculate  these  intersection  lengths, 
we  used  an  efficient  algorithm,  due  to  Siddon  [30],  in  which  the  intersections  of 
rays  with  separate  grids  of  parallel  horizontal  and  vertical  lines  are  computed, 
the  intersections  ordered  appropriately,  and  the  differences  computed.  In  the 
interest  of  storage  efficiency,  we  stored  the  weights  corresponding  only  to  a  single 
projection  view  at  any  one  time,  recomputing  them  anew  at  each  iteration.  To 
further  reduce  storage  requirements,  we  only  stored  non-zero  weight  elements 
along  with  a  separate  list  of  their  corresponding  indices. 

We  tested  the  algorithm  by  reconstructing  images  of  a  numerical  breast 
phantom  containing  large,  low-contrast  background  structures,  as  well  as  two 
small  lesions,  one  centrally  and  the  other  peripherally  located.  Because  all 
of  the  structures  are  ellipses,  the  phantom’s  projections  could  be  computed 
analytically,  and  we  generated  sinograms  of  64  angles  x  64  bins  and  16  angles  x 
64  bins.  We  reconstructed  each  of  these  using  ART  and  FBP  (using  a  ramp  filter 
with  cutoff  at  the  Nyquist  frequency).  We  then  added  Poisson  noise  (200,000 
total  counts  in  the  64-angle  sinogram,  50,000  counts  in  the  16-angle  sinogram) 
and  reconstructed  using  the  same  techniques. 

2.2.3  Results  and  conclusions 

The  results  of  the  ART  and  FBP  reconstructions  are  illustrated  in  Fig.  2. 


Figure  2:  Reconstructions  of  a  numerical  breast  phantom  using  ART  and  FBP. 
The  top  row  corresponds  to  a  64-angle  sinogram  and  and  the  bottom  row  to  a 
16-angle  sinogram.  Column  (a)  is  reconstructed  by  ART  from  noise- free  data, 
column  (b)  by  FBP  from  noise-free  data,  column  (c)  by  ART  from  noisy  data, 
and  column  (d)  by  FBP  from  noisy  data. 

The  ART  reconstructions  are  obviously  disappointing  and  especially  so  in  the 
presence  of  noise.  Both  methods  perform  well  for  the  noise-free  reconstructions 
from  64  angles,  but  in  the  presence  of  noise,  it  is  more  difficult  to  discern 
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the  lesions  in  the  ART  reconstruction  than  in  the  FBP  image.  Moreover,  the 
ART  reconstruction  contains  noise  structures  that  could  be  mistaken  for  lesions. 
In  the  noise-free,  16-angle  reconstructions,  the  lesions  are  again  more  visible 
in  the  FBP  reconstruction  despite  the  presence  of  distracting  star  artifacts. 
Finally,  while  the  lesions  are  difficult  to  discern  in  either  of  the  noisy  16-angle 
reconstructions,  the  FBP  reconstruction  again  seems  superior.  Coupled  with 
the  fact  that  ART  is  much  more  computationally  intensive  than  FBP,  these 
results  are  quite  disappointing. 

The  reasons  for  ART’s  poor  performance  in  the  noise-free  case  may  lie  in 
the  combination  of  numerical  instability  and  the  lack  of  regularization  in  the 
straightforward  ART  implementation  of  equation  (4).  Clearly  these  factors  are 
only  exacerbated  in  the  presence  of  noise.  These  poor  results  dissuaded  us  from 
pursuing  goal  3,  the  design  of  a  binary  ART  algorithm,  until  the  performance 
of  the  basic  iterative  algorithm  could  be  improved.  The  shortcomings  of  the 
ART  reconstructions  suggested  looking  into  means  of  regularizing  ART,  which 
in  turn  pointed  toward  ways  of  incorporating  explicit  statistical  information  into 
the  reconstruction  process,  a  step  that  had  been  planned  for  year  2  of  the  grant. 
At  present,  we  are  conducting  an  extensive  literature  search,  to  gain  exposure 
to  and  understanding  of  the  very  rich  field  of  statistical  image  reconstruction. 

2.3  Sinogram  interpolation  and  interpolation  method  prop¬ 
erties 

2.3.1  Background 

The  surprising  success  of  FBP  in  the  comparisons  with  ART,  even  for  few- view 
reconstructions,  prompted  us  to  reconsider  its  application  to  the  few- view  re¬ 
construction  problem.  It  is  well-known  that  the  principal  drawback  of  few- view 
FBP  reconstructions  is  the  presence  of  a  prominent  star-shaped  artifact  in  the 
background.  This  artifact  can  arise  even  when  the  angular  sampling  is  in  accor¬ 
dance  with  the  Nyquist  condition  and  is  due  to  the  fact  that  FBP,  being  based 
on  a  continuous,  analytic  Radon  inversion  formula,  implicitly  assumes  the  pro¬ 
jection  data  is  approximately  continuous  in  the  angular  variable.  If  the  data 
is  only  sparsely  sampled  in  the  angular  variable,  the  negative  side  lobes  intro¬ 
duced  by  the  reconstruction  filter  fail  to  cancel  in  the  intended  way,  giving  rise 
to  star-shaped  artifacts.  In  principal,  when  the  data  satisfies  the  Nyquist  condi¬ 
tion  at  least  approximately  but  fails  to  satisfy  FBP’s  assumptions — a  situation 
that  can  arise  when  the  object  being  imaged  is  relatively  simple  or  symmetric — 
additional  views  can  simply  be  interpolated  from  the  available  ones  and  FBP 
used  to  reconstruct. 

The  interpolation  method  used  should  of  course  exploit  the  periodicity  of 
the  data  in  the  angular  coordinate,  and  we  have  explored  four  popular  pe¬ 
riodic  interpolation  approaches:  linear  interpolation  with  periodic  boundary 
conditions,  periodic  spline  interpolation,  zero-padding  (ZP)  interpolation,  and 
circular-sampling  theorem  interpolation  (CST).  Linear  interpolation  needs  no 
introduction,  and  ZP  is  also  quite  familiar:  it  involves  extending  the  discrete 
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Fourier  transform  (DFT)  of  a  finite  sequence  with  zeroes  and  then  taking  an 
inverse  DFT  to  obtain  a  more  densely  sampled  version  of  the  sequence,  with  val¬ 
ues  interpolated  at  intermediate  positions  between  the  original  measured  sam¬ 
ples  [31-33].  CST  interpolation  is  a  special  case  of  Whittaker-Shannon  (W-S) 
sine  interpolation  that  applies  to  periodic  functions.  The  CST  holds  that  if  a 
function  g{x)  is  periodic  with  period  X,  it  can  be  reconstructed  from  a  finite 
number  N  of  samples  taken  over  one  period.  Specifically,  consider  g{x)  to  be 
sampled  such  that  gn  =  g{nAx),  where  Ax  =  X/N  and  n  =  0, . . . ,  iV  -  1.  If 
g{x)  is  bandlimited  to  frequency  K  (i.e.,  the  coeflScients  of  expansion  ak  of  the 
function’s  Fourier  series  satisfy  =  0  for  |A:|  >  K),  and  if  >  2K  -t- 1  and  is 
odd,  then  g{x)  can  be  reconstructed  exactly  from  its  samples  by  use  of 


AT-l 

9i^)  =  H  5n 

n=0 


{x  -  nAx)] 
AT  sin  [y  {x  —  nAx)]  * 


(5) 


Similarly,  if  N  >  2K  and  is  even,  then  the  value  of  g{x)  may  be  determined 
exactly  at  arbitrary  x  using 


q(x)-Tq  I  ^coc[^  (x  nAx) 

“is  1  ^sin[f  (x-nAx)] 


(6) 


As  for  periodic  spline  interpolation,  it  involves  fitting  the  data  with  pieceiwse 
cubic  polynomials  that  are  continuous  up  to  and  including  the  second  derivative 
at  the  joints  between  pieces  [34,35].  A  spline  g{x)  can  be  represented  by 

p(x)  =  a„  +  6„(x-x„)  +  c„(x-x„)^  +  d„(x-x„)®,  x  6  [x„,x„+i],  (7) 


and  it  can  be  shown  that  the  coefficients  On,  bn,  Cn,  and  dn  can  be  obtained 
from  the  vector  g  of  measured  data  points  through  matrix  multiplications  of 
the  form  a  =  Ag,  where  the  matrices  such  as  A  are  independent  of  the  data 
and  are  constructed  in  accordance  with  periodic  boundary  conditions. 

The  first  aspect  of  the  interpolation  work  was  to  explore  the  little-understood 
connections  among  ZP,  CST,  and  W-S  interpolation.  We  showed  that  they  are  in 
fact  mathematically  equivalent  for  the  task  of  increasing  the  density  of  angular 
samples,  with  zero-padding  being  considerably  more  computationally  efficient, 
so  we  could  omit  CST  from  the  list  of  approaches  being  considered  for  sino¬ 
gram  interpolation.  Next  we  used  statistical  hypothesis  testing  to  examine  the 
empirical  accuracy  of  the  interpolation  methods  for  the  task  of  interpolating  ad¬ 
ditional  sinogram  views.  Finally,  we  examined  the  noise  properties  of  the  various 
interpolation  methods,  deriving  expressions  for  the  variance  and  covariance  at 
points  in  interpolated  functions  given  that  the  measured  points  are  corrupted 
by  either  white  or  Poisson  noise. 
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2.3.2  Methods 

Connection  between  zero-padding,  CST,  and  W-S  interpolation 

The  connection  between  zero-padding,  CST,  and  W-S  interpolation  was  derived 
using  both  a  graphical  and  a  mathematical  approach  [36];  the  derivations  are 
different  for  the  cases  when  the  number  of  initial  samples  N  is  even  or  odd. 
In  the  interest  of  space,  we  summarize  here  only  the  graphical  proof  for  the 
case  when  N  is  odd.  Consider  an  arbitrary  complex  function  g{x),  which  is 
not  necessarily  periodic,  bandlimited,  or  compact,  and  which  has  a  Fourier 
transform  G{f).  The  real  parts  of  these  functions  are  illustrated  in  Fig.  3(a). 


Figure  3:  (a)  The  real  parts  of  an  arbitrary  function  g{x)  and  its  Fourier  trans¬ 
form  G{f).  Also  shown  is  the  subinterval  X  that  will  be  sampled,  (b)  The  DFT 
pair  samples  are  implicitly  samples  of  periodic  functions  in  both  the  spatial  and 
frequency  domains.  The  continuous  spatial- domain  function  g{x)  is  the  periodic 
extension  of  the  subinterval  X,  while  the  frequency-domain  function  G{f)  arises 
from  aliased  replications  of  the  original  Fourier  transform  further  modified  as 
a  result  of  the  truncation  of  g{x).  Shown  are  the  sampled  versions  of  g{x)  and 
G{f),  denoted  5s(x)  and  Gs(/),  respectively.  The  principal  periods  of  N  sam¬ 
ples  in  each  space  (which  correspond  to  the  N  values  of  the  DFT  pair)  are  also 
indicated,  (c)  The  operation  of  zero-padding  the  DFT  samples  is  equivalent  to 
multiplying  by  a  sequence  of  rect  functions,  as  shown. 

The  function  is  sampled  N  times  (where  N  is  an  odd  number)  over  a  sub-interval 
of  length  X,  with  sampling  distance  Ax  =  X/N  (without  loss  of  generality,  we 
assume  that  the  first  sample  lies  at  a:  =  0).  In  order  to  appreciate  the  effects 
of  zero-padding  the  DFT  of  this  sequence,  it  is  necessary  to  understand  the 
relationship  of  the  DFT  pair  to  the  true  function  g{x)  and  the  true  spectrum 
G{f).  The  DFT  pair  consists  of  the  N  space-domain  samples  of  g{x)  and  the 
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N  corresponding  frequency-domain  samples  yielded  by  the  DFT.  Implicitly, 
however,  each  sequence  represents  samples  of  a  single  period  of  a  continuous 
periodic  function  [37].  In  the  spatial  domain,  this  function,  which  we  call  g{x), 
is  simply  the  periodic  extension  of  the  subinterval  of  g{x)  that  we  have  sampled; 
that  is,  g{x)  is  the  periodic  function  that  arises  when  this  subinterval  of  length 
X  is  replicated  every  X  units.  In  the  frequency  domain,  we  have  a  function 
that  we  call  G(/),  which  is  periodic  with  period  1/Ax.  It  is  composed  of  aliases 
of  the  true  spectrum  and  is  generally  also  degraded  by  oscillations  introduced 
in  truncating  g{x).  The  sampled  functions,  which  we  denote  gsix)  and  G«(/), 
are  illustrated  in  Fig.  3(b)  and  can  be  represented  mathematically  as  gs{x)  = 

m  E"=-oo  and  Gsif)  =  G{f)  S{f^rn/NAx). 

Appreciating  the  implicit  periodicity  of  the  DFT  pair  is  important  in  un¬ 
derstanding  the  space-domain  effect  of  zero-padding  in  frequency  space  because 
any  modification  to  the  single  period  of  samples  available  affects  all  the  other 
periods  as  well.  With  this  in  mind,  it  is  easy  to  see  the  zero-padding  operation 
of  extending  the  N  available  frequency-space  samples  with  zeroes  to  PN  sam¬ 
ples  (where  P  is  an  integer)  is  equivalent  to  multiplying  the  sampled  periodic 
sequence  Gg{f)  by  a  periodic  sequence  of  rect  functions,  as  illustrated  in  Fig. 
3(c)  (P  =  2  in  this  illustration).  This  sequence  of  rects,  which  we  call  Z{f),  can 
be  represented  mathematically  as 


Z{f)  =  rect  (/Ax)  * 


(8) 


where  *  denotes  the  operation  of  convolution  and  where 

r  1  |/Ax|<  1/2 

rect(/Ax)  =  <  1/2  j/Ax|  =  1/2  (9) 

(  0  |/Ax|>  1/2. 

Thus  we  have  a  new  frequency-domain  function  G'(/)  =  Gs(/)Z(/),  and  by 
the  convolution-multiplication  theorem,  the  corresponding  new  spatial-domain 
function  g'{x)  is  given  by  ^'(x)  =  gs{x)  *z{x),  where  z{x)  is  the  inverse  Fourier 
transform  of  Z{f)  and  we  have  used  the  fact  that  ^s(x)  is  the  inverse  Fourier 
transform  of  G»(/).  Of  course,  z{x)  can  be  determined  from  Z{f)  by  another 
application  of  the  convolution-multiplication  theorem,  using  the  fact  that  the 
Fourier  transform  of  a  rect  is  a  sine  function  and  that  the  Fourier  transform  of 
a  sequence  of  delta  functions  is  itself  a  sequence  of  delta  functions.  Thus, 


g'i^)  = 


g{x)  S{x-nAx) 


sinc(a;/Ax)  ^  6  (x  - 

l=-oo  ^  ' 

It  can  be  shown  that  this  expression  is  equivalent  to: 

{oo  1  ^  /  Ax\ 

52  5(nAa;)sinc[(a: -nAa:)/Aa;]  V  52 

n=-oo  J  i=-oo  ^  ^ 


(10) 


(11) 
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which  indicates  that  zero-padding  the  DFT  of  the  Ax-spaced  samples  of  g{x) 
with  ZP  factor  P  is  equivalent  to  convolving  the  samples  of  g{x)  with  a  sine 
function  and  then  resampling  the  resulting  curve  with  spacing  AxfP. 

The  periodicity  of  ^(x)  makes  it  possible  to  collapse  the  infinite  sum  in  square 
brackets  to  a  sum  over  a  single  period  of  ^(x)  using  an  identity  that  lies  at  the 
core  of  the  CST.  Specifically, 


N-l 


Y:  5(nAx)smc  [(.  -  nAx)/Axl  =  Y  “if. _"„a1  j] 


n=0 


which  allows  us  to  write  equation  (11)  as 


(12) 


,'(x)  =  (x  -  if)  .  (13) 


The  proof  of  equation  (12)  is  given  by  Stark  [38]  in  the  context  of  his  derivation 
of  the  CST  (the  ratio  of  sine  functions  in  equation  (12)  is  the  CST  interpolation 
kernel),  and  while  the  functions  in  that  work  are  assumed  to  be  periodic  and 
bandlimited,  the  proof  of  the  above  identity  requires  only  the  assumption  of 
periodicity.  The  bandlimited  constraint  enters  only  in  the  claim  that  the  sine 
interpolation  on  the  left-hand  side  of  equation  (12)  is  exact,  a  claim  we  do 
not  yet  make.  At  this  point,  equations  (11)  and  (13)  simply  allow  us  to  claim 
that  using  zero-padding  with  zero-padding  factor  P  to  interpolate  from  an  odd 
number  N  of  Ax -spaced  samples  of  an  arbitrary  function  g{x)  is  equivalent  to 
applying  the  CST  interpolation  formula  to  those  samples  and  resampling  the 
resulting  curve  with  spacing  AxfP  and  is  also  equivalent  to  applying  the  W-S 
interpolation  formula  to  the  periodic  extension  of  those  samples  and  resampling 
that  resulting  curve  with  spacing  AxfP.  Cavicchi  [39]  derived  a  result  similar  to 
equation  (13)  without  explicitly  making  the  connection  to  the  CST.  Of  course, 
a  practical  implementation  of  zero-padding  only  yields  the  NP  samples  lying 
between  0  and  X.  If  ^(x)  is  actually  periodic  with  period  equal  to  the  length  of 
the  sampled  subinterval,  X,  then  ^(x)  =  ^(x),  and  if  g{x)  is  also  bandlimited 
to  the  Nyquist  frequency  1/2 Ax,  then  the  applications  of  the  W-S  and  CST 
formulas  in  equations  (11)  and  (13)  are  exact,  and  the  expression  in  brackets  in 
each  of  those  equations  is  simply  equal  to  ^(x).  Under  these  conditions,  then, 
both  equations  imply 


OO  y  A  \ 

g'{x)=9{x)  Y,  (14) 

/=-00  '  ' 

That  is,  under  these  conditions,  zero-padding  returns  exact  values  of  p(x)  spaced 
by  AxfP.  An  essentially  identical  result  was  derived  for  the  case  when  N  is 
even. 
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Emprical  accuracy  studies 

Having  established  this  connection  between  zero-padding  and  CST  interpola¬ 
tion,  and  their  exactness  for  periodic,  bandlimited  data,  we  then  examined  their 
accuracy  as  well  as  the  accuracy  of  periodic  spline  and  linear  interpolation  in  the 
face  of  typical,  noise-free  tomographic  data.  While  theoretical  arguments  can 
be  made  about  the  strengths  and  weaknesses  of  these  interpolation  approaches, 
the  arguments  generally  rely  on  assumptions  that  may  not  hold  exactly  for  typ¬ 
ical  tomographic  data.  In  view  of  this,  we  designed  an  empirical  comparison 
in  which  we  generated  100  different  sinograms  corresponding  to  a  family  of  nu¬ 
merical  breast  phantoms.  The  100  representatives  of  the  family  were  generated 
by  assigning  a  probability  distribution  to  the  parameters  specifying  the  ellipses 
that  made  up  the  phantom,  thus  resulting  in  a  varied  but  still  recognizable  set 
of  phantoms.  Noise-free  sinograms  of  512  angles  by  64  bins  were  generated  us¬ 
ing  the  analytical  expressions  for  the  Radon  transform  of  an  ellipse.  Each  of 
these  sinograms  was  subsampled  at  16,  32,  64,  128,  and  256  angles  and  then 
each  of  the  three  interpolation  approaches  was  used  to  interpolate  back  to  512 
angles.  The  normalized  root-mean-square  error  (NRMSE)  between  the  original 
and  interpolated  sinograms  was  computed.  Each  of  the  interpolated  sinograms 
was  then  reconstructed  using  FBP  and  the  NRMSE  between  the  reconstruction 
and  the  numerical  phantom  was  computed  as  well.  Finally,  paired  t-tests  and 
nonparametric  signed  rank  tests  were  used  to  assess  the  significance  of  any  dif¬ 
ference  in  the  performance  of  the  three  interpolation  approaches  from  each  of 
the  numbers  of  starting  samples. 

Noise  properties 

In  practice,  of  course,  the  samples  of  the  functions  of  interest  are  invariably  con¬ 
taminated  by  noise,  and  it  would  be  valuable  to  know  how  the  noise  corrupting 
the  measured  samples  is  propagated  into  the  interpolated  samples.  Consider  a 
periodic  function  g{x)  having  period  X  that  is  sampled  N  times  over  one  period, 
i.e.  at  points  Xn  =  Xn/N  for  n  =  0, . . . ,  -  1.  Assume  each  measurement  is 

corrupted  by  additive,  zero-mean  noise.  The  measured  samples  of  g{x)  can  then 
be  represented  as 

g(3:n)  =  (g(a:„))  +  n(x„),  (15) 

where  bold-faced  type  denotes  a  random  variable,  {)  represents  the  expectation 
operator  and  g(xn)  is  the  additive,  zero-mean  noise.  Let  g(a;)  be  the  curve  in¬ 
terpolated  from  the  noisy  samples  of  equation  (15).  We  compute  the  covariance 
of  this  function  for  two  points  x  and  x'  by  use  of 

cov(x,x')  =  ([g(x)  -  (g(x))]  [g(x')  -  (g(a:'))])  •  (16) 

Of  course,  once  the  covariance  has  been  computed,  the  variance  at  any  point 
x  is  given  by  var(x)  =  cov(x,rc).  Evaluation  of  equation  (16)  is  straightforward 
for  CST  and  ZP  interpolation,  but  more  subtle  for  spline  interpolation.  For  the 
first  two,  g(a:)  has  a  relatively  simple  form,  given  by 
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(17) 


where 


AT-l 

i(^)  =  ^  gix„)aN{x  -  Xn), 

n=0 


sin(iV6>/2)/7Vsin((9/2) 

sm[{N  -  l)e/2]/Nsin{e/2) 


Nodd 
N  even. 


(18) 


However,  for  periodic  splines,  as  we  see  in  equation  (7),  the  function  speci¬ 
fying  g(rc)  for  a  particular  x  depends  on  the  interval  in  which  x  falls.  This  can 
be  handled  by  defining  a  vector  of  functions,  indexed  by  the  intervals  between 
measured  samples  and  realizing  that  equation  (16)  will  yield  a  matrix  of  covari¬ 
ance  functions,  with  each  element  corresponding  to  a  pair  of  the  intervals  that 
X  and  x'  can  fall  into.  For  the  empirical  studies  we  calculated  analytically  128 
(and  127)  samples  of  a  typical  periodic  function  encountered  in  tomography:  the 
angular  function  corresponding  to  one  bin  of  the  projection  of  a  Shepp-Logan 
head  phantom.  We  generated  50,000  noise  realizations  of  these  samples  contam¬ 
inated  with  additive,  zero-mean  Gaussian  noise  and  50,000  others  with  Poisson 
noise.  We  then  interpolated  each  of  these  realizations  to  1024  (and  1016)  sam¬ 
ples  using  the  3  methods  discussed  above.  We  calculated  the  empirical  variance 
at  each  point  as  well  as  the  empirical  covariance  between  each  of  9  samples  and 
the  remaining  samples. 


2.3.3  Results  and  conclusions 

The  results  of  the  empirical  accuracy  study  can  be  summarized  as  follows.  In 
examining  the  sinogram  NRMSE  results,  the  following  ranking  held  for  both 
phantoms  and  for  virtually  every  number  of  starting  angles:  periodic  spline  in¬ 
terpolation  was  superior  to  linear  interpolation,  which  in  turn  was  superior  to 
zero-padding  interpolation.  The  same  general  trends  held  for  the  reconstructed 
image  NRMSE  results,  although  the  NRMSE  tended  to  change  very  little  for 
starting  numbers  of  angles  128  or  more,  the  errors  due  to  angular  sampling  and 
interpolation  effects  having  become  small  in  comparison  to  the  baseline  errors 
of  the  reconstruction  process  itself.  The  relatively  poor  performance  of  zero¬ 
padding  interpolation  comes  as  a  surprise,  especially  given  the  result,  reported 
above,  that  zero-padding  is  exact  and  equivalent  to  Whittaker-Shannon  inter¬ 
polation  when  applied  to  periodic,  bandlimited  data.  While  tomographic  data 
has  been  shown  to  be  at  least  quasi-bandlimited  in  the  angular  direction,  they 
are  not  exactly  bandlimited,  and  zero-padding  is  apparently  quite  sensitive  to 
departures  from  that  assumption.  Indeed,  on  inspecting  typical  zero-padding 
results,  we  saw  that  it  lead  to  ringing  artifacts  near  edges,  the  familiar  Gibbs 
phenomenon.  This  ringing,  which  extends  to  either  side  of  any  such  edge,  results 
in  a  relatively  large  NRMSE  compared  to  linear  and  spline  interpolation,  which 
react  more  locally.  In  conclusion,  then,  we  find  that  periodic  spline  interpola¬ 
tion  would  be  a  wise  choice  for  the  interpolation  of  additional  angular  views  for 


17 


noise- free  or  low  noise  data.  It  obviously  leads  to  smoother,  more  flexible  inter¬ 
polating  curves  than  linear  interpolation  with  fewer  of  the  non-local  responses 
to  violations  of  its  assumptions  that  plague  zero-padding. 

As  for  the  noise  studies,  we  derived  that  the  covariance  of  CST  and  ZP 
interpolated  functions  is 


N-lN-l 

cov(x,x')  =  XI  ”  Xj)(7N{x  -  Xk)  {n{xj)n{xk)) .  (19) 

j=0  k=0 

When  the  noise  is  white  (i.e.  (n(xj)ii(xfc))  =  crlSjk),  one  can  show  that 
cov(x,  x*)  =  (Tla^ix  -  x')*  This  has  the  interesting  consequence  for  the  variance 
that 


iVodd 
N  even. 


(20) 


We  note  two  interesting  facts  about  this  result.  First,  the  variance  is  constant 
at  all  points  in  the  interpolated  function,  including  points  corresponding  to  the 
original  measured  points.  Second,  in  the  case  when  N  is  even,  the  variance 
is  actually  reduced  relative  to  the  variance  in  the  measured  data  as  well  as 
to  the  variance  in  functions  interpolated  from  an  odd  number  of  points.  This 
argues  for  acquiring  an  even  number  of  samples  when  performing  this  kind  of 
interpolation,  though  the  difference  is  small  as  N  grows  large.  The  variance 
result  for  N  even  is  shown  in  Fig.  4(b). 

For  Poisson  noise  (i.e.  (n(xj)n(xA:)}  =  (p(iz^j))  one  can  show  that  for 
slowly  varying  functions  cov(x,x')  =  (g(x))  crp^{x  —  x')  and  thus  that  the  vari¬ 
ance  is 


(g(a^)) 

(e(*)>  (^) 


Wodd 

N  even. 


Again,  the  variance  remains  essentially  unchanged  from  that  found  in  the  sam¬ 
ples  except  for  the  slight,  uniform  reduction  when  N  is  even.  The  variance 
result  for  N  even  is  shown  in  Fig.  4(e). 

The  analytic  expression  for  the  spline  is  unwieldy  as  it  involves  a  sum  of 
many  products  of  the  matrices  like  A.  As  the  resulting  expression  must  be 
evaluated  numerically  to  discern  the  form  of  the  covariance  and  variance  curves, 
we  simply  present  the  result  of  this  evaluation  for  the  variance  for  N  even  and 
Gaussian  noise  in  Fig.  4(d),  along  with  empirical  variance  results  for  both 
Gaussian  and  Poisson  noise  in  Figs.  4(c)  and  4(f).  We  observe  here  that  the 
spline  interpolation  yields  reduced  variance  between  the  measured  samples,  with 
the  minimum  located  midway  between  samples  while  the  original  sample  points 
themselves  have  unchanged  variance.  N  being  odd  or  even  has  no  effect  on  the 
result. 

Overall,  then,  we  see  that  spline  interpolation,  in  addition  to  its  good  ac¬ 
curacy  also  has  desirable  variance  reduction  properties.  We  also  observe  that 
when  employing  CST  or  ZP  interpolation,  using  an  even  number  of  samples 
offers  a  slight  reduction  in  variance  at  no  cost  in  accuracy. 
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Figure  4:  (a)  Original  function  used  in  the  interpolation  studies.  All  the  results 
shown  here  involved  interpolating  from  128  to  1024  samples  of  this  analytically 
specified  function.  To  emphasize  detail,  we  show  the  variance  only  for  a  subin¬ 
terval  of  the  function,  (b)  Empirical  variance  for  CST  and  ZP  interpolation 
for  Gaussian  noise  ((7  =  4).  This  curve  matches  the  analytic  prediction  of  a 
constant  variance  equal  to  15.875.  (c)  Empirical  variance  for  periodic  spline  in¬ 
terpolation  for  Gaussian  noise  (cr  =  4).  The  maxima  correspond  to  points  that 
coincide  with  the  measured  samples,  (d)  Analytically  predicted  variance  for 
periodic  spline  interpolation  for  Gaussian  noise  (a  =  4).  (e)  Empirical  variance 
for  CST  and  ZP  interpolation  for  Poisson  noise.  The  variance,  while  locally  flat 
is  seen  to  track  the  variance  of  the  original  measured  samples,  (f)  Empirical 
variance  for  periodic  spline  interpolation  for  Poisson  noise.  The  variance  be¬ 
haves  locally  as  in  (c)  but  globally  is  seen  to  track  the  variance  of  the  measured 
samples. 

2.4  FBP  with  sinogram  preprocessing  and  effectively  2D 
smoothing 

2.4.1  Background 

Despite  the  mild  variance  reduction  offered  by  periodic  spline  interpolation, 
interpolating  among  noisy  data  is  still  a  risky  proposal.  In  principle,  it  would 
be  desirable  to  smooth  the  data  prior  to  interpolation.  Ideally,  we  would  like 
to  perform  a  fully  2D,  adaptive  smoothing  of  the  sinogram,  as  this  would  make 
maximal  use  of  the  correlations  inherent  in  the  data.  Adaptive  means  that 
the  degree  of  smoothing  applied  to  various  subsets  of  the  data  is  determined 
automatically  from  the  statistics  in  each  subset  separately,  in  contra^st  to  a 
global  smoothing  operation  that  applies  the  same  degree  of  smoothing  to  the 
entire  data  set.  However,  truly  2D  adaptive  smoothing  algorithms  are  in  general 
computationally  expensive  and  difficult  to  implement  [40].  In  response,  we  have 
developed  a  method  for  achieving  efficient,  effective  2D  adaptive  smoothing  by 
exploiting  Fourier  transforms  to  reduce  the  smoothing  problem  to  a  ID  one. 
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2.4.2  Methods 

Consider  a  2D  discrete  sinogram  p{^iyOj),  where  refers  to  the  ith  projection 
bin  {i  =  1, . . . ,  AT)  and  6j  the  jth  projection  angle  (j  =  1, . . . ,  M).  It  can  be 
shown  [41]  that  the  following  sequence  of  operations  is  equivalent  to  an  adaptive, 
2D  smoothing  of  the  sinogram: 

1.  Take  a  ID  discrete  Fourier  transform  of  the  sinogram  with  respect  to  the 
projection  angle  Oj;  the  result  can  be  viewed  as  a  set  of  ID  functions  of 
the  untransformed  variable  each  labeled  by  an  angular  frequency  index 
k,  i.e.  Pki^i). 

2.  Perform  an  adaptive  ID  smoothing  of  each  of  these  M  functions  of 

yielding  M  discrete  smoothed  functions  where  the  superscript  s 

stands  for  smoothed. 

3.  Perform  an  inverse  ID  discrete  Fourier  transform  of  P^i^i)  with  respect 
to  the  angular  frequency  k  to  recover  p^{^i,0j). 

The  adaptive  ID  smoothing  we  perform  on  each  of  the  M  functions  Pki^i) 
is  known  as  penalized  least-squares  smoothing  [27,28],  and  involves  fitting  the 
discrete  data  with  a  continuous  smoothing  curve  P^iO  that  minimizes  the  func¬ 
tional  ^ 

s{p^{0}  =  Y,[Pk{^i)  -  Pki^i)? + « rimorfd^,  (21) 

where  ^  is  a  continuous  variable  representing  the  position  along  a  given  projec¬ 
tion,  T  is  the  total  length  of  the  projection,  and  the  double  prime  denotes  the 
second  derivative  with  respect  to  The  two  terms  in  this  functional  represent 
the  competing  goals  of  achieving  a  good  fit  to  the  data  while  maintaining  a 
smooth  curve,  with  the  parameter  a  mediating  the  tradeoff.  For  instance,  if 
a  is  zero,  the  smoothness  constraint  disappears  and  the  minimizing  curve  will 
be  a  piecewise  linear  interpolant  to  the  data;  if  a  grows  large,  the  smoothness 
constraint  dominates  and  the  curve  approaches  a  simple  linear  fit  to  the  data. 
For  intermediate  values  of  a,  the  minimizing  curves  balance  the  goodness-of-fit 
and  smoothness  constraints.  It  can  in  fact  be  shown  that  the  minimizers  of 
this  functional  will  always  be  members  of  a  class  of  functions  known  as  natural 
cubic  splines.  These  are  piecewise  cubic  curves  that  join  at  the  abscissa  values 
where  they  are  continuous  up  to  and  including  the  second  derivative. 
Clearly  the  choice  of  a  determines  the  degree  of  smoothing  that  is  applied  to 
the  data,  and  it  is  this  parameter  that  is  determined  from  the  statistics  of  the 
data  itself  in  an  adaptive  implementation  of  penalized  least-squares  smoothing 
using  an  algorithm  known  as  generalized  cross-validation  [29].  Thus  a  generally 
different  a  is  used  in  the  smoothing  of  each  of  the  M  functions  Pk{^i)-  The 
resulting  continuous  smoothed  functions  PkiO  i^ust  then  be  sampled  to  yield 
the  discrete  functions  P^i^i)  which  are  subjected  to  the  inverse  DFT  in  step  3 
above. 
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It  should  also  be  noted  that  while  we  have,  for  the  sake  of  simplicity,  dis¬ 
cussed  effective  2D  smoothing,  the  technique  can  be  extended  to  any  number 
of  dimensions.  To  smooth  an  n-dimensional  function,  one  can  simply  take  an 
(n  -  l)-dimensional  Fourier  transform  of  the  function  and  then  perform  a  set 
of  ID  smoothings  over  the  untransformed  variable  prior  to  taking  an  inverse 
(n  “  l)-dimensional  Fourier  transform  [41], 

In  order  to  test  the  hypothesis  that  this  effectively  two-dimensional  adaptive 
smoothing  approach  is  superior  to  non-adaptive,  one-dimensional  smoothing 
such  as  the  apodization  used  in  filtered  backprojection,  we  again  calculated 
ideal-observer  SNRs,  this  time  for  images  reconstructed  by  FBP  after  various 
kinds  of  smoothing  were  applied  to  the  conventional  and  dedicated  SPECT 
projection  data  acquired  for  the  geometry  comparisons  described  in  Section  2.1. 
As  mentioned  there,  the  aim  here  is  to  see  which  smoothing  approach  best 
preserves  the  SNR  of  the  raw  projection  data.  The  detailed  procedure  is  as 
follows: 

1.  The  signal  projections  were  scaled  to  simulate  a  desired  tumor-background 
concentration  ratio  (6:1  in  this  case)  and  added  to  each  of  the  20  sets  of 
background  projections. 

2.  The  slice  through  the  center  of  the  lesion  was  selected  and  the  20  corre¬ 
sponding  signal-plus-background  sinograms  were  reconstructed  by  filtered 
backprojection  using  ramp  and  Hanning  filters  with  various  cutoff  frequen¬ 
cies  (0,4,  0.6,  0.8,  and  1.0  times  the  Nyquist  frequency).  The  sinograms 
were  also  processed  using  the  effective  2D  smoothing  procedure  described 
above  and  reconstructed  by  filtered  backprojection. 

3.  The  20  corresponding  sinograms  of  background  alone  were  processed  in 
the  same  way. 

4.  An  average  signal  image  was  determined  by  subtracting  the  average  of  the 
20  background  alone  reconstructions  from  the  average  of  the  20  signal- 
plus-background  reconstructions.  The  signal  power  spectrum  was  com¬ 
puted  by  squaring  the  Fourier  transform  of  this  image. 

5.  While  SPECT  images  are  not  stationary  in  general,  the  attenuated  pro¬ 
jections  of  a  uniform  cylinder  of  this  diameter  are  quite  flat  over  a  broad 
central  region,  and  thus  one  might  reasonably  expect  the  reconstructed 
images  of  this  cylinder  to  be  locally  stationary  near  their  center,  precisely 
where  the  lesion  is  expected  to  lie  [42].  The  “local”  Wiener  spectrum  in  this 
region  was  computed  from  the  20  images  of  background  alone  by  subtract¬ 
ing  the  average  background  image  from  each  of  the  individual  background 
images,  resulting  in  20  noise  images.  Each  such  image  was  multiplied  by 
a  circularly  symmetric  window  of  the  form: 

{1  lr|  <  0.9R 

I  {1  +  cos  [tt  (jrl  -  0.9i?)  /0.2R]}  0.9R  <  |rl  <  I.IR  (22) 
0  |r|  >  I.IR, 
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(with  appropriate  shifting  for  the  off-center  cylinder  in  the  conventional 
geometry),  where  R  is  the  radius  of  the  circular  region  over  which  the 
noise  is  expected  to  be  stationary  (chosen  to  be  6  pixels)  and  r  the  radial 
position  in  the  image.  The  power  spectrum  of  each  of  the  20  images  was 
computed  by  taking  the  square  of  the  Fourier  transform  of  the  resulting 
image.  The  20  power  spectra  were  then  averaged  and  scaled  so  that  the 
volume  under  the  Wiener  spectrum  equaled  the  average  variance  in  the 
circle  of  radius  R  [43-45]. 

6.  The  ideal  observer  SNR  was  then  determined  by  summing  the  quotient  of 
the  calculated  signal  spectrum  and  Wiener  spectrum. 

It  should  be  noted  that  in  using  the  ideal-observer  framework  at  all  it  is  implic¬ 
itly  being  assumed  that  the  data  satisfy  the  assumptions  discussed  in  Section 
2.1:  that  the  system  is  linear  and  that  the  noise  in  the  planar  or  reconstructed 
images  is  additive,  Gaussian,  zero-mean,  stationary,  and  independent  of  the 
presence  or  absence  of  the  signal.  Given  the  reasonably  high  count  levels  ("10- 
15 /pixel)  the  fact  that  the  signal  is  relatively  small  and  low  contrast,  and  the 
discussion  of  stationarity  in  point  5  above,  the  assumptions  about  the  noise 
seem  reasonable.  The  requirement  of  linearity  seemingly  undermines  the  use  of 
the  framework  to  analyze  images  that  have  been  processed  by  adaptive,  effec¬ 
tive  multi-dimensional  smoothing.  However,  what  is  truly  required  for  equation 
(1)  to  be  meaningful  is  not  linearity  in  the  face  of  any  possible  input  but  more 
specifically  that  the  system  transfer  function  be  the  same  whether  the  particular 
signal  of  interest  is  present  or  absent  from  the  particular  background  of  inter¬ 
est.  Again,  because  the  signal  in  question  is  relatively  small  and  low  contrast, 
it  should  not  greatly  affect  the  noise  properties  of  the  projection  images  and 
thus  the  effective  multi-dimensional  smoothing  algorithm  should  yield  a  similar 
effective  system  transfer  function  whether  or  not  the  signal  is  present. 

Having  verified  the  performance  of  the  smoothing  when  applied  to  sinograms 
containing  a  standard  number  of  projection  angles,  we  then  tested  it  on  few- 
view  sinograms,  after  which  periodic  spline  interpolation  was  used  to  generate 
additional  views  and  FBP  used  to  reconstruct.  We  acquired  projections  of  a 
ventricular  phantom  filled  with  3.27  mCi  of  Tc-99m.  The  phantom  was  not 
placed  within  a  water-filled  torso  phantom.  We  imaged  this  phantom  with  a 
Picker  PRISM  3000XP  three-headed  SPECT  system  fit  with  LEHRP  collima¬ 
tors.  We  acquired  SPECT  studies  containing  120  angular  views  over  360°  of 
the  phantom  placed  at  five  different  radial  offsets  from  the  center  of  rotation:  0, 
6,  9,  12,  and  15  cm.  We  used  a  25-cm  radius  circular  orbit  and  step-and-shoot 
mode  for  all  of  the  acquisitions;  each  head  acquired  to  a  128x128  pixel  image. 
A  total  of  about  500,000  counts  was  collected.  Prom  this  data  we  extracted 
sinograms  corresponding  to  a  single  transverse  slice  and  containing  15,  30,  60, 
and  120  views  respectively.  Thus  we  had  20  different  sinograms,  corresponding 
to  the  20  possible  combinations  of  radial  offset  and  number  of  angular  views. 
We  reconstructed  images  from  these  20  sinograms  in  3  ways:  (1)  Simple  recon¬ 
struction  from  the  available  views  by  conventional  FBP  using  a  Hanning  filter 
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(cutoff— 0.8).  (2)  Effective  2D  adaptive  smoothing  of  the  sinogram  and  recon¬ 
struction  from  the  available  views  by  FBP  using  a  ramp  filter  (cutoff=:=1.0).  (3) 
Effective  2D  adaptive  smoothing  of  the  sinogram,  followed  by  spline  interpola¬ 
tion  from  the  available  views  to  120  views,  and  reconstruction  by  FBP  using  a 
ramp  filter  (cutoffs  1.0). 

2.4.3  Results  and  conclusions 

The  ideal-observer  SNRs  for  the  detection  of  a  lesion  with  a  6:1  lesion-background 
concentration  ratio  are  listed  in  Table  2  for  the  two  different  tomographic  ge¬ 
ometries  and  various  kinds  of  smoothing.  For  reference,  we  also  include  the 
sinogram  SNR  values  discussed  in  Section  2.1  above.  The  results  indicate  that 


Smoothing  Method 

Conventional  SPECT 

Dedicated  SPECT 

Sinogram 

6.7 

10.7 

Hanning  (cutoff-0.4) 

2.4 

6.2 

Hanning  (cutoff-0.6) 

3.1 

8.0 

Hanning  (cutofF=0.8) 

3.5 

9.4 

Hanning  (cutoff— 1.0) 

3.7 

9.9 

Ramp  (cutoff=0,4) 

3.3 

8.1 

Ramp  (cutoff=0.6) 

3.7 

9.7 

Ramp  (cutoif=0.8) 

3.6 

9.9 

Ramp  (cutoff=1.0) 

3.7 

9.3 

Effectively  2D  smoothing 

5.5 

10.6 

Table  2:  Ideal-observer  SNRs  for  reconstructed  conventional  and  dedicated 
SPECT  images  of  a  6:1,  1-cm  lesion  using  various  kinds  of  smoothing. 


effective  two-dimensional  smoothing  provides  improvement  in  SNR  over  other 
filtering  approaches.  Images  corresponding  to  the  two  different  geometries  for 
selected  processing  methods  are  shown  in  Fig.  5. 

For  ease  of  comparison  and  simplicity  of  presentation  of  the  cardiac  phantom 
studies,  we  have  grouped  the  reconstructed  images  by  the  number  of  angles  in 
the  original  sinogram  and  we  show  in  Fig.  6  the  results  for  only  three  of  the 
radial  offsets:  0,  9,  and  15  cm.  We  observe  that  reconstructions  from  avail¬ 
able  views  without  pre-smoothing  or  interpolation  display  star-shaped  artifacts 
and  a  mottled  appearance  when  the  number  of  views  is  small.  2D  smoothing 
alone  reduces  the  noise  visibility  but  has  little  effect  on  the  star-shaped  arti¬ 
facts.  Effectively  2D  smoothing  followed  by  interpolation  to  120  views  reduces 
both  the  appearance  of  noise  and  the  star-shaped  artifacts.  However,  when  the 
radial  offset  is  large  and  the  number  of  initial  views  small,  interpolation  leads 
to  prominent  circular  artifacts.  By  considering  angular  sampling  requirements 
for  objects  of  corresponding  radial  extent,  it  can  be  shown  that  these  artifacts 
arise  precisely  when  the  initial  angular  sampling  is  insufficient.  No  processing 
can  compensate  for  this  sort  of  shortcoming  in  the  data.  However,  when  the 
sampling  is  adequate,  the  proposed  technique  allows  FBP  to  produce  high  qual- 
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Figure  5:  Reconstructed  slices  of  a  cylindrical  phantom  containing  a  1-cm,  6:1 
lesion  for  two  different  tomographic  geometries.  The  top  row  corresponds  to  the 
dedicated  SPECT  geometry,  the  bottom  row  to  conventional  SPECT.  The  left 
column  involves  FBP  reconstructions  with  a  Hanning  filter  (cutoff=0.8),  the 
middle  column  FBP  reconstructions  with  a  ramp  filter  (cutoff=0.6),  and  the 
right  column  FBP  reconstructions  after  effectively  2D  smoothing. 


ity  images  from  a  small  number  of  views.  Indeed,  we  have  conducted  further 
studies  with  a  cardiac  phantom  containing  a  lesion  defect,  and  found  that  for 
a  well-centered  phantom,  the  proposed  technique  yields  bullseye  plots  from  as 
few  as  15  or  30  angles  that  capture  the  defect  as  well  as  those  obtained  from 
120  views. 

2.5  Direct  spline  Radon  inverse 

2.5.1  Background 

Most  of  the  sinogram  preprocessing  approaches  described  so  far  have  involved 
fitting  of  continuous  smoothing  and  interpolating  splines  that  were  then  resam¬ 
pled  for  input  to  the  discrete  FBP  algorithm.  There  is  an  alternative  approach, 
however,  first  proposed  by  Wahba  [25]  that  allows  reconstruction  to  proceed 
directly  from  the  coefficients  of  the  splines  fit  at  each  projection  angle.  We  have 
simplified  her  approach,  using  relationships  among  spline  coefficients  to  elimi¬ 
nate  numerically  unstable  terms,  and,  perhaps  most  significantly,  extended  the 
approach  to  the  three-dimensional  Radon  inversion  problem,  where  it  takes  on 
a  much  simpler  form  than  in  two  dimensions  [46]. 

In  the  following,  we  will  denote  a  distribution  of  activity  in  two  dimensions 
by  f{x,y)  and  label  each  projection  through  it  by  the  pair  where  (f) 

specifies  the  projection  angle  and  ^  the  projection  distance.  The  value  of  such  a 
projection  is  given  by  the  line  integral  of  the  distribution  along  the  line  specified 
by  and  will  be  denoted  by  Similarly,  we  will  denote  a  three- 

dimensional  distribution  by  /(x,2/,z)  and  label  each  planar  projection  through 
it  by  the  triplet  {^,  0,  (^},  where  6  and  (p  specify  the  orientation  of  the  plane  and 
^  the  distance  of  the  plane  to  the  origin  of  the  coordinates.  The  value  of  such 
a  projection  is  given  by  the  planar  integral  of  the  distribution  over  the  plane 
specified  by  {^,  0,  (}>}  and  will  be  denoted  by  p(^,  0,  (p).  The  functions  p(^,  (p)  and 
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Figure  6:  Images  of  a  slice  through  a  Data  Spectrum  ventricular  phantom  recon¬ 
structed  by  FBP  after  the  indicated  processing  of  the  sinograms.  The  smoothing 
entailed  the  adaptive,  effectively  2D  smoothing  described  in  the  text,  while  in¬ 
terpolation  was  spline-based  interpolation  from  the  number  of  views  listed  to 
120  views. 


p{^,0,(j))  are  known  as  sinograms  because  in  the  two-dimensional  case  a  point 
distribution  in  {rr,  y]  space  maps  to  a  sinusoid  in  {^,  (f)}  space. 

The  Radon  transform  is  a  continuous,  integral  transform  that  relates  the 
sinogram  to  f{x^y)  in  two  dimensions  and  the  sinogram  p{^,0,(j))  to 

f{x,y^z)  in  three  dimensions  [29,47,48].  Inverting  the  Radon  transform  ex¬ 
actly  to  recover  a  distribution  requires  continuous,  noise-free  knowledge  of  the 
distribution’s  sinogram,  which  entails  having  an  infinite  set  of  perfect  projec¬ 
tion  measurements.  In  practice,  of  course,  one  can  only  collect  projection  data 
in  the  two-dimensional  case  for  a  finite  number  of  projection  distances  (we 
call  these  discrete  samples  projection  bins)  at  a  finite  number  of  projection  an¬ 
gles  (l)j,  and  the  measurements  are  invariably  contaminated  with  noise.  In  the 
three-dimensional  case,  the  planar-integral  projection  data  cannot  generally  be 
measured  directly  and  must  instead  be  generated  from  line-integral  projection 
data;  it  is,  however,  still  only  generated  for  a  finite  number  of  projection  bins 
and  projection  angles  (j)j  and  Ok- 

Because  the  sinogram  p(^,  (p)  in  two  dimensions  (or  p(^,  0,  (p)  in  three  di¬ 
mensions)  is  known  only  on  a  finite  grid  of  and  (pj  (or  (pj^  and  Ok)^  we 
cannot  invert  the  Radon  transform  exactly  to  recover  the  distribution  /(x,y) 
(or  f{x,y,  z))\  we  must  instead  turn  to  a  discrete  approximation  of  the  inverse. 
For  instance,  one  way  of  implementing  the  most  popular  two-dimensional  Radon 
inversion  algorithm — filtered  backprojection  (FBP)  [49, 50] — begins  with  a  dis- 
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Crete  filtration  of  the  sinogram.  The  filtered  samples  of  for  each  projection 
angle  are  then  backprojected  onto  the  image  grid  and  the  resulting  images 
summed  to  give  a  final  reconstructed  image.  One  way  to  view  the  backprojec- 
tion  step  is  to  imagine  casting  a  perpendicular  ray  from  each  image  pixel  to  each 
projection  angle  in  turn,  summing  the  sinogram  values  picked  up  at  each  angle 
to  obtain  the  final  pixel  value.  In  this  view,  the  difficulty  lies  in  determining 
what  value  to  pick  up  from  each  projection  angle,  for  in  general  the  perpendic¬ 
ular  line  will  not  fall  directly  in  the  center  of  a  projection  bin.  In  the  simplest 
schemes,  one  simply  picks  up  the  value  of  the  bin  the  pixel  projects  onto,  while 
in  a  more  complicated  approach  one  might  perform  a  linear  interpolation  of  the 
values  in  the  two  nearest  projection  bins.  In  the  most  sophisticated  schemes, 
one  uses  a  weighted  average  of  the  values  in  a  slightly  larger  neighborhood.  A 
similar  procedure  can  be  used  to  implement  three-dimensional  FBP  [51]. 

The  discreteness  of  the  sinogram  thus  dictates  discreteness  in  both  the  fil¬ 
tration  and  backprojection  steps  of  the  algorithm.  If,  however,  one  had  a  con¬ 
tinuous,  analytic  expression  for  the  sinogram  at  each  projection  angle — if  the 
sinogram  were  a  set  of  known  one- dimensional  functions  of  ^ — it  might  be  pos¬ 
sible  to  implement  the  filtration  and  backprojection  steps  in  a  continuous  man¬ 
ner.  The  filtration  could  be  performed  analytically,  and  the  resulting  filtered 
projections  would  be  continuous  functions  which,  in  the  pixel-traversal  view  of 
backprojection  described  above,  could  be  sampled  wherever  a  projection  might 
strike  without  any  need  to  interpolate.  Naturally,  such  an  analytic,  continuous 
expression  for  the  sinogram  cannot  be  obtained  directly  from  any  tomographic 
imaging  system,  but  must  rather  be  obtained  by  fitting  an  analytic  expression 
to  the  discretely  sampled  data.  Not  every  class  of  function  that  could  be  fit  to 
the  data  would  allow  the  filtration  of  the  projections  to  be  calculated  in  closed 
form,  but  Wahba  [25]  has  shown  that  one  class  that  does  allow  such  a  solution 
are  the  cubic  splines,  piecewise  third-order  polynomials  that  are  continuous  up 
to  and  including  the  second  derivative  at  the  connection  points  between  pieces. 
This  is  the  class  of  fitting  functions  we  investigate  in  this  paper,  introducing 
Wahba’s  results  (with  some  corrections  and  simplifications),  and  extending  the 
treatment  to  the  three-dimensional  Radon  transform.  This  method  offers  a  cer¬ 
tain  conceptual  appeal,  as  well  as  the  advantage  of  directness  when  one  wishes 
to  smooth  noisy  projection  data  by  fitting  curves  that  minimize  the  popular 
penalized  least-squares  measure  [35].  As  it  turns  out,  the  minimizers  of  this 
measure  are  the  natural  cubic  splines  mentioned  above  and  thus  reconstruction 
can  proceed  directly  from  the  coefficients  of  these  splines  in  this  case. 

2.5.2  Methods 

Inverse  2D  Radon  transform  in  coordinate  space 

The  essential  problem  in  two-dimensional  tomography  is  the  reconstruction  of 
a  distribution  f{x,y)  (or  /(r,  0)  in  polar  coordinates)  from  knowledge  of  the 
discrete  sinogram  p(^i,0j),  where  i  =  and  j  =  1,...,M.  This 

convention  for  the  index  z,  particularly  the  choice  of  an  odd  number  of  projection 
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bins,  will  simplify  the  mathematical  expressions  to  be  derived  below,  but  the 
proposed  technique  is  applicable  to  geometries  with  an  even  number  of  bins  ajs 
well  We  assume  in  the  present  paper  that  we  have  a  parallel-beam  geometry. 

Perhaps  the  most  familiar  way  of  expressing  the  inverse  of  the  two-dimensional 
Radon  transform  is  in  terms  of  the  frequency-space  representation  of  the  con¬ 
tinuous  sinogram: 

f{r,0)=  /  \v\P{v,<l>)e^^”'^^dvd<i>,  (23) 

where  v  is  the  spatial-frequency  variable  corresponding  to  P{v,(j>)  is  the 
Fourier  transform  of  the  sinogram  with  respect  to  the  variable  and  j  is  the 
imaginary  number  This  expression  provides  the  theoretical  basis  for  FBP, 
as  lt;|  is  just  the  familiar  ramp  filter.  This  expression  may  be  written  in  coordi¬ 
nate  space  as 


where 


fir,9)  =  ^I^JrA4>)d4>, 

=  jinn  + 


(24) 

(25) 


and  in  which  =  r  cos(^  -  <j))  and  p'(^,  0)  is  the  first  derivative  of  the  sinogram 
with  respect  to  (  [47,52].  Taking  the  limit  as  e  0  of  the  sum  of  these 
two  integrals  allows  us  to  avoid  integrating  over  the  singularity  at  In 

general  equations  (24)  and  (25)  are  less  useful  than  equation  (23)  because  the 
data  collected  in  PET,  SPECT,  or  CT  constitute  samples  of  the  sinogram  itself, 
and  do  not  provide  any  direct  information  about  p'(C,  0).  However,  by  fitting 
an  analytic,  differentiable  function  of  ^  to  the  projection  data  at  each  angle,  we 
could  obtain  an  expression  for  If  had  an  auspicious  functional 

form,  we  would  then  be  able  to  solve  the  integrals  in  equation  (25)  in  closed 
form. 


Interpolating  and  smoothing  splines 

To  obtain  an  expression  for  the  sinogram  that  is  analytic  and  differentiable  with 
respect  to  the  variable  we  fit  a  function  to  the  projection  data  at  each  angle. 
That  is,  for  each  angle  0j,  we  fit  a  one-dimensional  function  of  p^,.  (0,  to  the 
sinogram  values  (^i)  measured  on  the  2N-\-l  abscissas  {i  =  -‘N, . . . ,  N). 
If  the  data  is  noiseless,  it  is  desirable  to  use  a  function  that  passes  through  the 
points  P0j  (^i),  which  we  call  an  interpolating  curve,  while  if  the  data  is  noisy 
a  smoothing  curve  may  be  more  appropriate.  One  fitting  framework  that  can 
handle  both  of  these  situations  is  known  as  penalized  least-squares  [35],  in  which 
the  function  p^p.  {^)  is  chosen  to  be  the  minimizer  of  the  functional 

=  E  b^,(^i)-P0,(Ci)]'+A (26) 

i=-N  “ 
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where  a  and  b  are  the  endpoints  of  the  interval  on  which  the  curve  (^)  is  to  be 
defined.  The  first  term  in  equation  (26)  is  the  familiar  squared-error  measure, 
while  the  second  is  a  measure  of  the  smoothness  of  the  fit  curve.  The  parameter 
A  thus  mediates  the  tradeoff  between  the  competing  goals  of  achieving  a  good 
fit  to  the  data  and  maintaining  a  smooth  curve.  By  choosing  A  to  be  zero, 
we  eliminate  the  smoothness  constraint  and  ensure  that  the  minimizing  curve 
will  be  an  interpolant  to  the  data;  if  A  grows  large,  the  smoothness  constraint 
dominates  and  the  curve  approaches  a  linear  fit  to  the  data.  For  intermediate 
values  of  A  the  minimizing  curve  balances  the  goodness-of-fit  and  smoothness 
constraints.  The  parameter  A  can  be  chosen  a  priori  [34]  or  it  can  be  determined 
from  the  statistics  of  the  data  using  an  automatic  procedure  such  as  generalized 
cross-validation  (GCV)  [53|. 

The  minimizers  of  this  functional  F  belong  to  the  class  of  functions  known 
as  natural  cubic  splines  [34].  These  are  piecewise  polynomial  curves  that  join  at 
the  abscissa  values  where  they  are  continuous  up  to  and  including  the  second 
derivative.  They  can  be  represented  as 


P4>s  (0  =  +  Cie/2  +  dif/S,  C  €  [^i,  6+i] ,  (27) 

where  a^,  6^,  Ci,  and  di  are  constants  that  fully  specify  the  spline  curve  on  the 
interval  [Cz,^i+i].  Of  interest  to  the  Radon  inversion  problem  is  the  fact  that 
we  can  approximate  the  first  derivative  of  the  sinogram  for  fixed  angle  (/>j  by 

p'(e,  <t>i)  =f{0  =  bi  +  Ci^  +  die,  4  e  [6,  ii+i]  ■  (28) 

Spline-based  inverse  of  the  2D  Radon  transform 

Given  this  analytic  expression  for  the  derivative  of  the  sinogram,  we  can  proceed 
with  the  inversion  of  the  two-dimensional  Radon  transform  in  equation  (24). 
While  the  sinogram  now  has  a  continuous  representation  in  the  variable  it 
is  still  discrete  in  the  angular  variable.  Assuming  that  the  M  angular  samples 
are  equally  spaced  over  180°  or  360°  (the  result  is  the  same  in  either  case),  the 
integral  in  equation  (24)  can  be  approximated  by  the  sum 

1  ^ 

=  (29) 

j=l 

where  Jr,e{<t>j)  is  given  by  equation  (25).  For  a  given  coordinate  (r,0)  in  image 
space,  and  for  a  given  projection  angle  0j,  =  rcos(0  —  0^).  We  label  the 
projection  bin  that  falls  in  by  m,  that  is,  G  [^m,  Cm+i]-  Using  the  expression 
for  p' {(,(!>)  given  by  equation  (28),  Jr,e{<j>j)  can  be  expressed  as 


N-l 

E 

i=— TV,  i^m— 


hi 


+ Cj^ + die 
i'-i 
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(30) 


+  lim 

C-+0 


If 


^  ,  |*^m+i  _j_  ^  dm^ 


-d^  + 


rK^+ 


d^ 


where  the  integrals  of  equation  (25)  are  now  expressed  as  sums  of  integrals 
over  the  subintervals  between  the  original  abscissa  points,  with  appropriate 
accommodation  for  the  singularity  at  These  integrals  can  be  solved  in  closed 
form  [25],  and  the  resulting  expression  contains  some  potentially  unstable  terms. 
However,  by  combining  these  terms  in  a  particular  way  and  invoking  spline 
identities,  a  stable  form  can  be  derived: 

JrA<t>j)  =T+  ibi  +  Cii'  +  di^)  In  (/_~/*  ) 

N~l  .  N-\ 

-  (c,  +  2d,e')(6+i  -  ^i)  +  2  E  J  ’ 


i=-N 


i=-N 


where  T  is  given  by 

T  =  (6„_1  +  Cm-lC  +  In  ie  -  U-l) 

~  (^>m+l  +  Cm+1^^  +  dm+1^'^)  ln(^m+2  ~  ^0-  (32) 


The  3D  Radon  transform  in  coordinate  space 

The  essential  problem  in  three-dimensional  computed  tomography  is  the  re¬ 
construction  of  a  distribution  f{x,y,z)  from  knowledge  of  the  discrete  planar- 
integral  sinogram  0j),  where  i  =  —N,  =  1, . . . ,  and  k  = 

1, . , . ,  M6i  [54].  In  general,  these  planar  integrals  are  not  measured  directly  by 
tomographic  imaging  systems,  but  must  rather  be  calculated  by  “rebinning^’  the 
line  integrals  that  are  measured  directly  [51].  The  inverse  three-dimensional 
Radon  transform  has  a  form  similar  to  the  two-dimensional  case,  with  a  few  dif¬ 
ferences  that  greatly  simplify  the  task  of  evaluating  it  numerically.  Specifically, 

fix,  y,z)  =  -^  £  p"  4>)  sin  ed4>de,  (33) 


where 


^'  =  x  sin  6c,os(j)-\-y  sin  z  cos 


(34) 


and  p”  (^',  0, 0)  is  the  second  derivative  of  the  three-dimensional  sinogram  with 
respect  to  evaluated  at  ^  This  expression  differs  in  two  principal  ways 
from  the  expression  for  the  two-dimensional  inverse  Radon  transform  given  by 
equations  (24)  and  (25).  First,  it  now  involves  the  second  derivative  of  the 
sinogram  with  respect  to  ^  rather  than  the  first  derivative  and  second,  the 
convolution  in  ^  has  disappeared.  This  refiects  the  fact  that  an  inverse  Radon 
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transform  of  odd  degree  can  be  calculated  using  purely  local  information — the 
value  of  the  image  at  a  point  can  be  determined  solely  from  information  at  points 
in  the  sinogram  space  that  (x,  y,  z)  projects  onto,  rather  than  from  a  convolution 
integral  over  all  points  in  sinogram  space  as  in  the  even-dimensional  case  [48,55]. 

Spline-based  inverse  of  the  3D  Radon  transform 

As  in  the  two-dimensional  case,  we  wish  to  fit  an  analytic  function  of  ^  to  each 
sequence  of  labeled  by  a  distinct  pair  {(j)j,0k}-  We  do  so  using  the  spline 
formalism  described  above  and  obtain 


P4>j,Ok{0  —  /S,  ^  *  (3^) 

Consequently  the  second  derivative  of  the  sinogram  is  approximately 

p"  =  Ci  +  2di^,  ^  .  (36) 

Now  as  in  the  two-dimensional  case,  the  discreteness  of  the  angular  samples 
means  that  the  two  integrals  give  way  to  sums  and  we  write 

Me 

/(x, y,  z)  =  — sin0/.,  (37) 

^  j=l  k—l 

where 

Jx,y,z{0k,4>j)  =  Ci  +  2di^',  ^6[Ci,^i+l],  (38) 

and 


=  X  sin  Ok  cos  0  j  H-  y  sin  Ok  sin  <f)j  +  z  cos  Ok . 

The  simplicity  of  the  three-dimensional  inversion  is  now  apparent,  for  the 
functions  J  can  be  evaluated  in  a  straightforward  manner  whereas  in  the  two- 
dimensional  case  evaluation  of  the  functions  J  involved  performing  a  compli¬ 
cated  integral  over  ^  and  taking  care  in  handling  numerically  unstable  terms. 

Application  to  phantom  and  real  data 

In  order  to  demonstrate  that  the  2D  direct-spline  inverse  of  the  Radon  transform 
produces  images  of  comparable  or  superior  quality  to  filtered  backprojection 
(FBP),  we  reconstructed  images  of  a  numerical  Hoffman  brain  phantom  [56] 
using  both  methods.  The  sinogram  consisted  of  128  simulated  noiseless  pro¬ 
jections  of  the  phantom,  taken  over  360°  and  each  comprising  400  projection 
bins.  The  sinogram  contained  a  total  of  1.72x10®  counts.  We  first  reconstructed 
the  phantom  using  standard  area- weighted  FBP  with  a  ramp  filter  (cutoff=1.0 
times  the  Nyquist  frequency).  We  then  fit  an  interpolating  spline  to  the  data 
at  each  of  the  128  projection  angles  and  used  the  2D  direct-spline  technique 
to  reconstruct.  Poisson  noise  was  then  added  to  the  sinogram.  Images  were 
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reconstructed  using  FBP  with  a  ramp  filter  (cutofF=1.0)  as  well  as  the  direct- 
spline  technique  using  an  interpolating  spline  in  order  to  see  how  the  algorithms 
compared  in  the  presence  of  noise  without  prior  smoothing.  Smoothing  splines 
were  then  fit  to  the  projection  data  at  each  angle,  using  GCV  to  determine  the 
smoothing  parameter,  and  an  image  reconstructed  from  the  coefficients  using 
the  direct-spline  method.  In  order  to  examine  the  performance  of  FBP  in  the 
face  of  data  with  the  same  degree  of  smoothness,  we  sampled  these  smoothing 
splines  to  generate  a  discrete  sinogram  that  was  reconstructed  using  FBP  with 
a  ramp  filter  (cutoff  =1.0). 

For  the  3D  case,  we  reconstructed  images  of  a  Data  Spectrum  ventricu¬ 
lar  phantom  from  projection  data  acquired  on  a  Picker  XP3000  three-headed 
SPECT  system  fitted  with  high-resolution,  parallel-hole  collimators.  The  phan¬ 
tom  was  filled  with  3.27  mCi  of  Tc-99m  and  placed  at  the  center  of  rotation. 
Each  head  collected  data  on  a  128x128  grid  and  at  120  projection  angles  over 
360°.  We  rebinned  the  projection  data  from  a  single  head  to  generate  planar 
integrals  on  a  128x60x120  grid.  Images  were  reconstructed  from  this  planar- 
integral  data  using  3D  FBP,  and  also  using  the  3D  direct-spline  inversion  method 
after  splines  were  fit  to  the  data  as  described  above.  We  then  fit  smoothing 
splines  to  the  planar-integral  data  and  reconstructed  directly  from  the  spline 
coefficients.  Finally,  we  sampled  the  smoothing  splines  to  generate  a  smoothed, 
discrete  sinogram  and  used  that  as  input  to  the  3D  FBP  algorithm. 

Resolution,  noise,  and  signal-to-noise  studies 

In  order  to  compare  quantitatively  the  resolutions  of  the  direct-spline  algorithms 
(both  two-  and  three-dimensional)  with  their  FBP  counterparts,  we  acquired 
projection  images  of  a  small  (1-cm)  spherical  lesion  containing  7.6  mCi  of  Tc-99m 
placed  in  an  800  cc  cylindrical  phantom  containing  cold  (zero-activity)  water. 
A  Picker  XP2000  two-headed  SPECT  system  fitted  with  ultra-high-resolution, 
parallel-hole  collimators  was  used.  The  heads  rotated  at  their  minimum  radius 
of  rotation  (9  cm)  and  acquired  120  views  over  360°  onto  a  128x128  matrix 
(pixel  size=4.67  mm). 

For  the  two-dimensional  algorithms,  we  extracted  the  2D  sinogram  corre¬ 
sponding  to  the  slice  through  the  center  of  the  lesion  and  reconstructed  images 
using  FBP  with  a  ramp  filter  (cutoff=1.0)  as  well  as  using  the  direct  spline 
inversion  with  interpolating  splines.  The  reconstructed  lesion  was  approxi¬ 
mately  a  symmetric  2D  Gaussian  in  shape  and  we  determined  its  full-width 
half-maximum  by  collapsing  it  into  a  one-dimensional  function  and  fitting  this 
profile  with  a  Gaussian  curve.  For  the  three-dimensional  reconstruction  algo¬ 
rithms,  we  rebinned  the  projection  data  to  generate  planar-integral  data  on  a 
128x60x120  grid.  We  then  used  3D  FBP  and  the  3D  direct  spline  method  (using 
interpolating  splines)  to  reconstruct  the  slice  through  the  center  of  the  lesion 
and  determined  the  FWHM  of  the  resulting  Gaussian  by  the  same  method  as 
above. 

In  order  to  isolate  the  contribution  of  the  reconstruction  algorithm  to  the 
FWHM  of  the  lesion  in  the  reconstructed  images,  the  contribution  from  the 
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lesion’s  inherent  width  as  well  as  the  imaging  system’s  point-spread  function  had 
to  be  removed.  The  net  effect  of  these  two  factors  was  estimated  by  determining 
the  average  FWHM  of  the  lesion  as  it  appeared  in  the  120  projection  images. 
Assuming  then  that  the  reconstruction  algorithms  could  be  characterized  by 
Gaussian  point-spread  functions,  the  FWHM  of  these  functions  were  determined 
by  subtracting  (in  quadrature)  the  average  projection  FWHM  from  the  FWHMs 
of  the  reconstructed  lesions  discussed  above. 

To  characterize  the  noise  level  in  images  reconstructed  by  the  direct-spline 
methods  and  their  FBP  counterparts,  we  acquired  20  1-minute  projection  data 
sets  of  the  same  800  cc  cylindrical  phantom  used  in  the  resolution  studies  above, 
this  time  containing  3.7  mCi  of  Tc-99m  and  no  lesion.  One  slice  of  this  uniform 
cylinder  was  reconstructed  for  each  of  the  20  datasets  using  2D  FBP,  3D  FBP, 
2D  spline  inversion,  and  3D  spline  inversion  (both  of  these  using  interpolating 
splines).  For  a  given  algorithm,  the  same  six  circular  regions  of  interest  (ROI) 
were  examined  in  each  of  the  20  slices  and  the  coefficient  of  variation  (the 
standard  deviation  of  the  pixel  values  in  the  ROI  divided  by  the  mean  of  the 
pixel  values  in  the  ROI)  calculated  for  each  of  the  120  ROIs.  The  average  of 
these  120  coefficients  of  variation  was  then  computed. 

It  is  not  uncommon  for  a  reconstruction  algorithm  to  offer  enhanced  reso¬ 
lution  at  the  price  of  amplified  noise.  The  overall  effect  of  such  a  tradeoff  is 
sometimes  better  characterized  by  computing  a  signal- to-noise  ratio  (SNR).  We 
used  the  two  datasets  described  above  to  compute  the  ideal-observer  SNR  using 
the  method  described  in  section  2.4.2. 

2.5.3  Results  and  conclusions 

The  results  of  reconstructing  the  Hoffman  brain  phantom  with  and  without 
noise  using  both  the  2D  direct-spline  inverse  and  2D  FBP  are  depicted  in  Fig. 
7.  The  algorithms  clearly  yield  qualitatively  similar  results.  The  results  of 
reconstructing  the  ventricular  phantom  data  using  3D  direct-spline  inversion 
with  interpolating  splines,  3D  FBP,  3D  direct-spline  inversion  with  smoothing 
splines,  and  3D  FBP  using  a  sinogram  resampled  from  the  smoothing  splines 
are  depicted  in  Fig.  8.  The  algorithms  are  again  seen  to  yield  qualitatively 
similar  results. 

The  resolution  measurements  for  the  four  basic  algorithms — 2D  direct-spline 
inversion,  2D  FBP  with  a  ramp  filter  (cutoff— 1.0),  3D  direct-spline  inversion, 
and  3D  FBP — are  summarized  in  Table  3.  The  results  indicate  that  the  direct- 


Algorithm 

FWHM 

2D  direct  spline 

1.6  mm 

2D  FBP 

4.5  mm 

3D  direct  spline 

3.9  mm 

3D  FBP 

5.0  mm 

Table  3:  FWHM  of  in-plane  reconstruction  point-spread  functions 
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Figure  7:  Reconstructions  of  a  Hoffman  brain  phantom  from  simulated  projec¬ 
tions  without  noise  (a-b)  and  with  Poisson  noise  added  (c-f).  Reconstruction 
methods  are:  (a)  FBP  with  ramp  filter  (cutoff— 1.0),  (b)  direct-spline  inversion 
using  interpolating  splines,  (c)  FBP  with  ramp  filter  (cutoff— 1.0),  (d)  direct- 
spline  inversion  using  interpolating  splines,  (e)  FBP  from  a  sinogram  gener¬ 
ated  by  sampling  smoothing  splines,  (f)  direct-spline  inversion  using  smoothing 
splines. 


spline  inversions  have  superior  resolution  to  FBP  in  both  the  2D  and  3D  cases 
and  also  that  the  2D  algorithms  have  superior  resolution  to  the  3D  algorithms. 
The  results  of  the  noise  study  are  summarized  in  Table  4,  where  it  can  be 


Algorithm 

COV 

2D  direct  spline 

0.60 

2D  FBP 

0.39 

3D  direct  spline 

0.34 

3D  FBP 

0.23 

Table  4:  Coefficients  of  variation  for  various  reconstruction  algorithms 

seen  that  the  noise  level  in  the  direct-spline  reconstructions  is  higher  than  that 
in  the  FBP  reconstructions  and  that  the  noise  level  in  the  2D  reconstructions 
is  higher  than  that  in  the  3D  reconstructions. 

Finally,  the  ideal-observer  signal-to-noise  ratio  results  are  summarized  in 
Table  5  for  the  four  basic  algorithms  as  well  as  their  counterparts  using  smooth¬ 
ing  splines.  We  observe  that  the  direct-spline  algorithms  have  slightly  lower 
SNRs  than  their  FBP  counterparts  when  using  interpolating  splines,  while  the 
SNRs  become  comparable  when  using  smoothing  splines.  Furthermore,  the  use 
of  smoothing  splines  seems  to  have  little  effect  on  SNR  in  the  2D  case  while 
degrading  it  in  the  3D  case.  All  of  the  reconstructed  images  are  seen  to  have 
lower  SNRs  than  the  raw  projection  data,  a  fact  that  is  discussed  in  greater 
detail  in  the  next  section.  Typical  images  reconstructed  using  each  of  these 
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Figure  8:  Selected  slices  of  a  ventricular  phantom  reconstructed  by  (a)  3D  direct- 
spline  inversion  using  interpolating  splines,  (b)  3D  direct-spline  inversion  using 
smoothing  splines,  (c)  3D  FBP,  and  (d)  3D  FBP  from  a  sinogram  obtained  by 
sampling  the  smoothing  splines  in  (b). 


Algorithm 

Ideal-observer  SNR 

Sinogram 

10.7 

2D  interpolating  spline 

7.9 

2D  FBP 

8.6 

2D  smoothing  spline 

8.4 

2D  FBP  w/  smth.  spline 

8.4 

3D  interpolating  spline 

9.0 

3D  FBP 

9.6 

3D  smoothing  spline 

7.5 

3D  FBP  w/  smth.  spline 

7.0 

Table  5:  Ideal-observer  SNRs  for  various  reconstruction  algorithms 


eight  methods  are  shown  in  Fig.  9. 

As  discussed  in  section  2.5.1,  the  principal  difference  between  the  2D  and  3D 
direct-spline  inversion  algorithms  and  FBP  is  in  the  nature  of  the  interpolation 
step  upon  backprojection.  The  interpolation  in  FBP  is  simply  less  accurate 
than  the  more  sophisticated  cubic-spline  interpolation  used  in  the  direct-spline 
method.  The  cruder  FBP  interpolation  is  more  likely  to  smooth  over  high- 
frequency  variations  in  the  projection  data  than  is  cubic-spline  interpolation, 
and  thus  it  is  not  surprising  that  the  FBP  algorithms  have  inferior  resolution 
to  the  spline-based  algorithms,  as  illustrated  in  Table  3.  However,  because  the 
high-frequency  components  of  the  data  include  considerable  noise  as  well,  the 
FBP  algorithms  would  be  expected  to  produce  less  noisy  reconstructions  than 
the  spline-based  reconstructions.  This  expectation  is  confirmed  by  the  results 
of  the  noise  study  reported  in  Table  4, 

The  3D  direct-spline  and  FBP  algorithms  are  both  seen  to  have  inferior  reso¬ 
lution  and  lower  noise  levels  than  their  2D  counterparts.  This  can  be  attributed 
to  the  fact  that  the  3D  reconstruction  process  involves  an  additional  averaging 
or  smoothing  step  which  occurs  when  the  raw  projection  data  is  rebinned  into 
a  planar-integral  sinogram  by  performing  area-weighted  forward  projections  of 
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Figure  9:  Reconstructions  of  a  selected  slice  of  a  cylindrical  phantom  containing 
a  spherical  lesion.  Reconstruction  methods:  (a)  2D  direct-spline  inversion  using 
interpolating  splines,  (b)  2D  direct-spline  inversion  using  smoothing  splines,  (c) 
2D  FBP,  (d)  2D  FBP  from  a  sinogram  obtained  by  sampling  the  smoothing 
splines  in  (b),  (e)  3D  direct-spline  inversion  using  interpolating  splines,  (f)  3D 
FBP,  (g)  3D  direct-spline  inversion  using  smoothing  splines,  (h)  3D  FBP  from 
a  sinogram  obtained  by  sampling  the  smoothing  splines  in  (g) 


the  2D  projection  data  at  each  projection  angle. 

Table  5  lists  the  ideal-observer  SNR  for  the  raw  projection  data  and  for  the 
various  reconstruction  approaches.  It  is  a  fact  that  processing  or  even  image 
reconstruction  can  never  improve  the  ideal-observer  SNR  over  that  found  in  the 
raw  projection  data.  However,  these  operations  can  certainly  diminish  the  SNR 
if  they  are  in  some  way  singular  and  if  the  signal  vector  has  components  in  the 
null  space.  We  observe  that  all  the  reconstructed  image  SNRs  are  in  fact  lower 
than  that  of  the  raw  projection  data.  The  differences  among  those  reconstructed 
image  SNRs  give  some  clue  as  to  how  much  of  the  information  contained  in  those 
projections  is  preserved  in  the  reconstructed  image.  For  instance,  we  see  that 
ideal-observer  SNRs  are  slightly  lower  for  the  2D  and  3D  direct-spline  inversions 
using  interpolating  splines  than  for  their  FBP  counterparts  using  ramp  filters. 
For  this  particular  detection  task,  then,  the  interpolating-spline  algorithm’s  am¬ 
plification  of  noise  outweighs  the  improvement  it  affords  in  resolution  relative  to 
FBP.  However,  when  the  noise  is  mitigated  prior  to  reconstruction,  as  when  the 
projection  data  has  been  fitted  with  smoothing  splines,  the  SNR  gap  between 
the  spline  algorithms  and  FBP  is  considerably  narrowed  in  the  2D  case  and 
reversed  in  the  3D  case.  This  suggests  that  the  spline-based  algorithms  may  be 
of  greatest  use  when  resolution  is  paramount  and  the  data  contains  relatively 
little  noise,  a  situation  more  often  encountered  in  computed  tomography  than 
in  nuclear  medicine.  The  use  of  smoothing  splines  is  seen  to  provide  little  or  no 
improvement  in  SNR  in  the  2D  case.  Smoothing  does  affect  the  ideal-observer 
SNRs  in  the  3D  case,  but  for  the  worse.  It  is  clear  from  examining  the  images 
of  Fig.  9  that  the  reconstructed  images  using  3D  algorithms  and  smoothing 
splines  have  an  oversmoothed  appearance.  This  can  most  likely  be  attributed 
to  the  fact,  mentioned  above,  that  the  3D  reconstruction  process  effectively  in¬ 
volves  a  prior  smoothing  during  the  rebinning  step.  The  adaptive  smoothing 
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algorithm  certainly  makes  allowances  for  the  lower  variability  in  the  rebinned 
data  and  smooths  this  data  less  than  it  would  the  raw  projection  data.  How¬ 
ever,  the  modified  statistics  of  the  rebinned  data  simply  do  not  agree  as  well 
with  the  statistical  model  assumed  by  the  smoothing  algorithm,  so  it  is  perhaps 
not  surprising  that  it  yields  a  sub-optimal  result.  It  remains  a  topic  for  fur¬ 
ther  investigation  as  to  whether  smoothing  the  raw  projection  data  prior  to  the 
rebinning  step  produces  a  better  result. 


3  Conclusions 

The  work  of  the  past  year  has  touched  on  a  wide  variety  of  topics  related  to  the 
reconstruction  of  dedicated  SPECT  scintimammographic  images  from  a  smaller 
number  of  views  than  is  usually  used  in  nuclear  medicine  tomography.  The 
following  principal  conclusions  have  emerged  from  the  work: 

1.  An  ideal-observer  SNR  analysis  of  the  data  obtained  in  planar,  conven¬ 
tional  SPECT,  and  dedicated  SPECT  scintimammography  indicated  that 
the  dedicated  SPECT  geometry  is  the  best  of  the  three  for  detecting  typ¬ 
ical  malignant  breast  lesions  given  equal  total  imaging  time  for  all  the 
modalities.  This  finding  reaffirms  the  assumption  underlying  this  project: 
that  despite  the  fact  that  no  one  has  yet  developed  a  dedicated  SPECT 
sctintimammographic  system,  its  obvious  advantages  make  it  almost  in¬ 
evitable  that  it  will  win  a  role  in  the  diagnosis  and  management  of  breast 
cancer. 

2.  Further  application  of  the  ideal-observer  framework  to  study  the  depen¬ 
dence  of  dedicated  SPECT  scintimammography  SNR  on  the  number  of 
projection  views  acquired  revealed  an  approximate  square  root  relation¬ 
ship.  This  implies,  for  instance,  that  cutting  in  half  the  number  of  pro¬ 
jection  views  acquired  only  results  in  a  v/2  reduction  in  the  SNR. 

3.  Images  of  a  numerical  breast  phantom  reconstructed  by  an  unregularized 
ART  algorithm  were  found  to  be  clearly  inferior  to  those  reconstructed  by 
ramp-filtered  FBP,  and  especially  so  in  the  presence  of  noise.  Moreover, 
even  in  the  case  of  a  few- view  reconstruction,  the  lesions  were  more  de¬ 
tectable  in  the  FBP  reconstruction,  despite  prominent  star  artifacts,  than 
they  were  in  the  ART  reconstruction.  This  finding  had  three  consequences 
for  the  research.  First,  it  dissuaded  us  from  pursuing  the  implementation 
of  a  binary  version  of  the  ART  algorithm  until  adequate  regularization 
approaches  have  been  explored.  Second,  it  impelled  us  to  begin  exploring 
such  regularization  approaches  as  well  as  the  related  statistical  image  re¬ 
construction  approaches.  Third,  it  inspired  us  to  reexamine  FBP  and  to 
explore  combinations  of  smoothing  and  interpolation  algorithms  for  sino¬ 
gram  preprocessing  in  the  hopes  of  eliminating  the  star  artifact  common 
in  few- view  FBP. 
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4.  The  first  conclusion  to  emerge  from  our  study  of  interpolation  methods 
was  the  mathematical  equivalence  of  zero-padding  and  circular-sampling 
theorem  interpolation  and  the  little-appreciated  fact  that  zero-padding  is 
exact  and  equivalent  to  Whittaker-Shannon  interpolation  when  interpo¬ 
lating  a  periodic,  bandlimited  function  from  a  set  of  samples  satisfying  the 
Nyquist  condition.  This  equivalance  allowed  us  to  omit  CST  interpolation 
from  the  set  of  those  being  considered  for  use  in  sinogram  interpolation, 
as  zero-padding  has  much  more  favorable  computational  properties. 

5.  Studies  of  the  noise  properties  of  interpolation  methods  yielded  three  in¬ 
teresting  conclusions.  The  first  is  that  when  interpolating  from  samples 
corrupted  by  white  noise,  CST  and  ZP  yield  interpolated  points  with  con¬ 
stant  variance.  The  second  is  that  when  the  number  N  of  initial  samples  is 
odd,  this  constant  variance  equals  the  variance  in  the  measured  samples, 
while  when  the  number  of  initial  samples  is  even  it  equals  (A^  - 1)  /N  times 
the  variance  in  the  measured  samples.  This  indicates  that  there  is  a  slight 
advantage  to  using  an  even  number  of  samples,  at  least  when  N  is  small. 
The  third  is  that  periodic  spline  interpolation  yields  nonconstant  variance 
in  the  interpolated  points:  equal  to  the  initial  variance  at  points  coincident 
with  measured  points  and  falling  to  a  minimum  between  such  points.  This 
suggests  that  periodic  spline  interpolation  has  a  slight  variance-reducing 
advantage  over  ZP  interpolation. 

6.  Studies  of  the  empirical  accuracy  of  various  types  of  periodic  interpolation 
approaches  using  statistical  hypothesis  testing  lead  to  the  conclusion  that 
periodic  spline  interpolation  is  superior  to  linear  or  zero-padding  interpo¬ 
lation  in  the  face  of  typical  tomographic  data.  This  finding,  coupled  with 
the  variance-reduction  property  discussed  in  point  5  indicate  that  peri¬ 
odic  spline  interpolation  is  the  best  choice  for  interpolation  of  additional 
angular  views  in  tomography. 

7.  We  developed  an  effectively  two-dimensional,  adaptive  smoothing  approach 
that  exploits  Fourier  transforms  to  reduce  the  dimensionality  of  the  smooth¬ 
ing  problem  and  we  used  ideal-observer  analysis  to  show  that  use  of  this 
approach  prior  to  reconstruction  by  FBP  preserves  more  of  the  informa¬ 
tion  content  of  the  sinogram  than  the  use  of  spatially  invariant  apodization 
filters.  The  technique  was  then  applied  to  few- view  data  and  then  followed 
by  interpolation  of  additional  views  with  promising  results. 

8.  Finally,  we  examined  and  tested  a  direct  spline-based  Radon  inversion 
technique  that  minimizes  the  inaccuracies  introduced  by  the  usual  dis¬ 
cretization  of  the  image  reconstruction  process.  While  computationally 
intensive  in  two  dimensions,  the  technique  is  quite  straightforward  in  the 
three-dimensional  case.  In  both  cases  it  was  found  to  yield  higher  res¬ 
olution,  albeit  noisier  images  than  FBP,  with  slightly  reduced  SNRs.  It 
remains  to  be  seen  whether  it  enhances  human  observer  performance. 
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